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I.  INTRODUCTION 


LSPRWC  is  a  FORTRAN  77  computer  program  written  to  perform  a  least  squares 
polynomial  regression  on  a  given  set  of 

(x,  y)  data  pairs  or  observations.  In  common  parlance,  LSPRWC,  a  “curve  fit”  procedure,  finds 
that  np  order  polynomial 


Pipe)  =  ^  bjxi 
j= o 

which  minimizes  the  sum  of  the  squared  errors,  the  squared  residuals,  between  the  polynomial 
value  P(x)  and y(x)  over  the  set  of  (x,  y)  data  pairs. 

In  addition,  LSPRWC  provides  for  the  imposition  of  constraints  on  P(x)  such  that  the  nc 
derivative  of  P(x)  may  be  specified  at  a  given  x  location. 
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LSPRWC  was  written  for  LINUX  /UNIX  users  as  an  analysis  tool  in  science  and 
engineering.  The  program  was  intended  for  users  who  care  about  numerical  precision,  wish  to 
have  a  well  documented  source  code  in  hand,  and  abhor  spreadsheet  calculations  and  license 
fees. 


LSPRWC  requires  command  line  interaction.  No  graphical  user  interfaces  are  provided 
with  the  computer  program — although  they  could  certainly  be  added.  The  code  was  not  meant 
for  the  “point  and  click”  community  of  computer  users. 

Least  squares  polynomial  regression  is  a  well  known  analytical  technique  dating  back  to  at 
least  the  18th  century3.  The  technique  remains  a  premier  method  for  the  functional 
approximation  of  presumably  repeatable  phenomena  when  observations  of  that  phenomena  are 
known  to  be  imperfect.  The  technique  is  also  commonly  used  to  quantify  those  phenomena 
which  either  defy  description  by  the  known  laws  of  physics  or  for  which  the  mathematical 
descriptions  are  simply  too  unwieldy  or  numerically  intensive  for  practical  application. 

Although  the  functional  form  for  linear  regression  need  not  be  based  on  polynomials, 
polynomials  have  proven  to  be  particularly  useful  for  a  broad  range  of  applications. 

Additionally,  the  closed-form  mathematical  solution  for  P{x)  is  readily  derived,  as  given  in 
Appendix  A.  Once  the  bj  coefficients  are  known,  polynomial  evaluation  reduces  to  a  simple 
numerical  process.  Furthermore,  it  is  relatively  easy  to  compute  derivatives  and  other  properties, 
such  as  arc  length  and  integrals  for  polynomials. 

LSPRWC  was  specifically  written  for  source  code  control,  source  code  that  is  readily 
tailored  for  specific  applications,  and  to  maintain  ease  of  use  with  the  advent  of  new 
computational  platfonns  and  operating  systems.  As  such,  the  code  eliminates  those  inexplicable 
fee  requirements  for  licensed  software  which  exploits  an  analytical  technique  embodied  in  the 
public  domain  for  centuries. 

The  addition  of  constraints  was  added  to  LSPRWC  as  a  particularly  useful  capability  not 
usually  included  with  commercial  least  squares  regression  software  packages.  Often  as  not, 
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some  physical  constraint — an  initial  condition  such  as  zero  velocity  and  a  second  derivative  such 
as  gravitational  acceleration — are  known  absolutely  and  can  be  included  with  LSPWRC  to 
better  approximate  a  given  phenomena. 

LSPRWC,  as  written  for  LINUX/UNIX  operating  systems,  requires  only  a 
FORTRAN  774  compiler  and  a  tenninal  window  system  such  as  the  X  Window  System5 .  The 
code  is  known  to  compile  successfully  with  the  Free  Software  Foundation  G776  compiler  and  to 
work  well  on  Apple  Mac  OS X1  and  Yellow  Dog8  LINUX  computer  operating  systems.  Although 
written  to  FORTRAN  77  standards,  the  code  may  compile  with  FORTRAN  959  compilers  and 
does  so  with  the  Intel  ifort10compiler.  LSPRWC  memory  and  disk  storage  requirements  are 
flexible  and  readily  adjusted. 

The  LSPRWC  computer  code  is  not  guaranteed  to  run  successfully  or  operate  without 
failure  on  any  given  computer  and/or  operating  system.  Considerable  error  messaging  is 
included  in  the  FORTRAN  77  code,  but  this  does  not  assure  that  the  linear  regression  process 
operates  as  intended.  Hence,  it  is  advised  that  the  user  first  execute  test  cases  with  known 
solutions  for  comparison  with  the  LSPRWC  code  output. 

II.  FEATURES 

The  LSPRWC  code  was  written  to  provide  a  user  friendly  environment  while  offering 
extensive  analytical  capabilities  by  means  of  the  following  features: 

•  Full  directives  are  written  to  the  terminal  window  whenever  input  to  the  LSPRWC 
code  is  requested. 

•  All  operational  parameters  are  input  using  either  NAMELIST  or  free  field  reads. 

•  Input  to  the  LSPRWC  code  from  formatted  disk  files  may  be  offered  as  an  option. 
In  such  cases,  the  data  file  fonnat,  including  free  field  reads,  will  be  input  from  the 
command  line. 

•  Sets  of  (x,  y )  data  pairs  are  entered  either  by  NAMELIST  or  read  from  a  formatted 
disk  file. 

•  Optional  constraints  are  input  either  by  NAMELIST  or  read  from  a  formatted  disk 
file.  Each  constraint  is  defined  by  the  derivative  order  nc  =  0  through  np-  1,  the  x 
location  for  the  constraint,  and  the  value  of  the  nc  derivative  of  P(x)  at  the  constraint 
location. 

•  Computational  results  from  the  linear  least  squares  regression  are  output  to  the 
tenninal  window  as  the  LSPRWC  code  executes.  This  output  includes  the  values 
of  the  bj  coefficients  of  the  np  order  polynomial  P(x)  and  the  standard  deviation  of 
the  ( x ,  y)  data  pairs  with  respect  to  the  regression  polynomial. 

•  As  an  option,  x  values  may  be  prescribed  for  evaluation  of  the  resultant  least 
squares  polynomial.  These  x  values  may  be  entered  by  NAMELIST,  input  from  a 
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formatted  disk  file,  or  alternatively  as  minimum,  maximum,  and  interval  values  to 
produce  an  array  of  evaluated  (x,  P(x))  data  pairs. 

•  A  plot  option  is  included  in  the  LSPRWC  code  such  that  the  user  may  visually 
compare  the  evaluated  (x,  P(pc))  data  with  the  (x,  y)  data  pairs.  No  plot  package  is 
actually  included  with  the  LSPRWC  code,  but  the  two  data  sets,  (x,  P(x))  and  (x,y), 
are  written  to  fonnatted  data  files  for  ready  access  to  the  user’s  favorite  plot 
package.  The  plot  option  subroutine  QPS  can  also  easily  be  modified  by  the  user 
to  include  an  in-line  plot  routine  for  execution  with  the  LSPRWC  code. 

•  Output  of  results  from  a  LSPRWC  run  may  optionally  be  written  to  a  fonnatted 
disk  file  as  either  the  entire  set  of  run  parameters,  input  data,  and  results  or  as 
simply  the  evaluated  data.  In  the  latter  case,  the  evaluated  data  includes  the  choices 
of  x,  the  0  through  np  -  1  derivative  value  of  P(x),  the  radius  of  curvature,  and  the 
arc  length  as  columnated  data  sets. 

•  Restart  options  are  included  in  the  LSPRWC  code  allowing  for  easy  modification 
of  operational  parameters,  for  example  polynomial  order  np,  without  having  to  re¬ 
enter  all  of  the  input  information. 

•  Constraints  are  satisfied  using  the  method  of  undetermined  Lagrange 
Multipliers11. 

•  The  LSPRWC  code  includes  a  matrix  inversion  subroutine  GJEMPS  which  is 
based  on  Gauss-Jordan  elimination  using  a  maximum  pivot  strategy  .  With  some 
minimal  effort,  the  user  may  substitute  their  preferred  matrix  inversion  technique 
should  that  be  desired. 

•  All  of  the  FORTRAN  77  array  sizes  are  set  in  PARAMETER  statements  within 
the  LSPRWC  code  making  it  relatively  easy  to  find  and  change  those  array  size 
limits  should  that  be  required. 

•  The  LSPRWC  code  maintains  64-bit  precision  with  IMPLICIT  REAL*8  statements 
in  each  subroutine  as  required. 

III.  SPECIAL  FEATURES 

The  following  special  features  were  added  to  the  LSPRWC  code  to  enhance  the  user 
friendly  environment: 

•  If  a  simple  Yes/No  input  is  requested,  then  any  upper/lower  case  variation  of  Yes, 
No,  Y,  or  N  provides  a  permissible  entry  as  controlled  by  the  FORTRAN 
subroutine  YNOUS  listed  in  Appendix  C. 
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Disk  files  to  be  used  for  LSPRWC  code  input  may  contain  alphanumeric  headers 
such  as 


Thermal  Data  Set 


T(K) 

289.26111111111 

289.81666666667 

290.37222222222 

290.92777777778 

291.48333333333 

292.03888888889 

292.59444444444 

293.15000000000 


P(Pa) 

0.18305580366000E+04 
0.18967477057200E+04 
0. 19650058020000E+04 
0.20346428497200E+04 
0.21070378003200E+04 
0.2181501 1780800E+04 
0.22587224587200E+04 
0.23373226908000E+04 


The  disk  tile  data  will  be  read,  line  by  line,  and  ignored  until  the  prescribed  data 
format  is  satisfied. 

•  Command  line  NAMELIST  input  can  be  a  nuisance.  NAMELIST  entries  begin  in 
column  2  with  either  a  $  or  &  sign  followed  immediately  with  the  NAMELIST 
name.  The  entries  continue  with  the  name  of  the  variables  to  be  changed  from  their 
current  values,  an  equal  sign,  and  the  new  value  to  be  assigned  to  the  variable.  The 
NAMELIST  entry  must  then  end  with  either  a  $  or  &  sign.  The  FORTRAN 
subroutines  NDRUS  and  NDWUS,  as  given  in  Appendix  C,  were  included  in  the 
LSPRWC  code  to  simplify  this  process. 

If  a  NAMELIST  input  is  requested,  then  all  variables  are  written  to  the  tenninal 
window  with  default  values  to  facilitate  cut  and  paste  input  using  a  three-button 
mouse. 

The  starting  column  for  NAMELIST  entries  is  irrelevant. 

NAMELIST  entries  need  not  begin  with  the  $  or  &  sign  nor  is  the  NAMELIST 
name  required. 

The  NAMELIST  entry  must  end  with  either  the  $  or  &  signs  but  only  if  variable 
values  are  altered.  A  simple  keyboard  return  will  accomplish  a  NAMELIST  entry 
with  no  changes  to  the  variable  default  values. 

•  Whenever  FORTRAN  statements  include  an  option  for  redirection  following  an 
error,  e.g. 


READ(UNIT=4,NML=PARM,ERR=169), 

and  the  source  of  the  error  is  trivial,  for  example  a  misspelled  NAMELIST  variable 
name,  then  provision  is  included  in  the  LSPRWC  code  for  correction  of  the  error 
without  the  need  to  restart. 
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IV.  PROGRAM  STRUCTURE 


The  structure  of  the  LSPRWC  FORTRAN  77  computer  code  as  listed  in  Appendix  B  is 
reasonably  straight  forward.  The  program  routines  and  their  purposes  are  as  follows: 

MAIN — Executive  control  for  program  LSPRWC 

XYDIS — Subroutine  for  input  of  the  (x,  y)  data  set  or  observations 

CDS — Subroutine  for  input  of  the  constraints 

PODS — Subroutine  for  input  of  the  polynomial  order 

LS  P  R  W  C — Subroutine  to  perform  the  least  squares  polynomial  regression  with 
constraints 

XDPEDS — Subroutine  to  provide  an  x-coordinate  set  for  evaluation  of  the  least  squares 
regression  polynomial 

PES — Polynomial  evaluation  subroutine 

QPS — Quick  plot  subroutine 

FDFOS — Formatted  disk  fde  output  subroutine 

GJEMPS — Matrix  inversion  subroutine  using  Gauss-Jordan  elimination  with  maximum 
pivot  strategy 

POWER — Subroutine  to  calculate  integer  powers  of  real  numbers 
FACTRL — Factorial  evaluation  subroutine 

As  can  be  seen,  the  actual  least  squares  polynomial  regression  algorithm  resides  in 
subroutine  LSPRWC.  Subroutines  CDS,  PODS,  XDPEDS,  and  XYDIS  provide  or  define 
input  infonnation  for  program  LSPRWC.  Subroutine  PES  provides  output  data  for  program 
LSPRWC.  Subroutines  FDFOS  and  QPS  output  information  from  program  LSPRWC. 
Subroutines  FACTRL,  GJEMPS,  and  POWER  are  computational  utility  routines.  A  flow  path 
for  program  LSPRWC  is  given  in  Figure  1. 
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Additional  low-level  utility  routines  are  listed  in  the  FORTRAN  Tool  Library  in 
Appendix  C. 

V.  LIMITATIONS 

The  LSPRWC  code  has  the  following  requirements  or  limitations: 

•  The  LSPRWC  code  was  written  with  the  following  set  of  COMMON  and 

PARAMETER  statements  which  set  the  array  sizes  consistently  in  all  subroutines. 
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MNXYDP  =  Maximum  Number  of  (x,  y )  Data  Pairs 
=  101 

MNXCE  =  Maximum  Number  of  x-Coordinates  for  polynomial  evaluation 
=  101 

MXANC  =  Maximum  Number  of  Allowable  Constraints 
=  20 

MXNR  =  Maximum  Number  of  Rows 
=  30 

MXNRM1  =  Maximum  Number  of  Rows  minus  1 
=  29 

MXNC  =  Maximum  Number  of  Columns 
=  31 

MXNCM1  =  Maximum  Number  of  Columns  minus  1 
=  30 

The  number  of  rows  and  columns  refers  to  the  size  of  the  augmented  matrix  for  the 
solution  of  n  linear  equations  in  n  unknowns.  It  follows  from  the  derivation  of 
Appendix  A  that 


MXNR  > 

rip  +  £  +  1 , 

(1) 

MXNRM1  = 

MXNR  -  1, 

(2) 

MXNC  = 

MXNR  +  1 ,  and 

(3) 

MXNCM1  = 

MXNC  -  1. 

(4) 

The  PARAMETER  variables  MXANC  and  MXNR  apply  only  for  NAMELIST 
input.  Any  number  of  ( x ,  y )  data  pairs  or  x-coordinates  for  polynomial  evaluation 
can  be  read  from  disk  files. 

•  The  set  of  ( x ,  y)  data  pairs,  as  read  into  program  LSPRWC,  must  be  monotonically 
increasing  in  x.  Provision  could  have  been  made  within  program  LSPRWC  to 
reorder  any  set  of  (x,  y)  data  pairs;  however,  it  has  been  discovered  that  the  failure 
of  a  data  set  to  satisfy  this  requirement  is  often  the  indication  of  unknown  errors. 
Hence,  no  data  set  reordering  has  been  implemented. 

•  The  number  of  (x,  y)  data  pairs  or  observations,  m,  must  satisfy 
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m  >  np  +  C  +  1  . 


(5) 


•  The  polynomial  order  np  must  satisfy 

0<np<9.  (6) 

The  LSPRWC  code  could,  though  not  readily,  be  modified  for  polynomial  orders 
greater  than  9;  however,  this  largely  defeats  the  purpose  intended  for  polynomial 
regression. 

•  The  constraint  order  nc  must  satisfy 

0  <  nc  <  np  .  /yx 

•  Polynomials  P(x)  will  not  produce  an  infinite  slope.  Hence,  a  transformation — 
possibly  a  rotation — may  be  required  for  a  given  set  of  (x,  y )  data  pairs  before  a 
least  squares  polynomial  regression  can  be  completed  to  eliminate  that  possibility. 

VI.  ERROR  MESSAGES 

The  LSPRWC  code  was  written  with  the  following  extensive  embedded  error  messaging: 

•  The  error  messages  cover  the  standard  FORTRAN  errors  and  those  associated 
directly  with  the  operation  of  program  LSPRWC. 

•  Wherever  possible,  the  LSPRWC  code  identifies  the  source  of  the  particular  error 
and  the  subroutine  call  structure  leading  to  that  error. 

•  Should  readily  correctable  errors  occur,  for  example  a  misspelled  NAMELIST 
variable  name  entry,  then  program  LSPRWC  will  identify  the  error  and  loop,  in 
this  case  for  the  NAMELIST  entry. 

•  The  standard  FORTRAN  errors  source  from  the  FORTRAN  OPEN,  CLOSE, 
READ,  and  WRITE  commands. 

•  Program  LSPRWC  will  issue  an  error  message  should  the  array  limits  be  exceeded. 
(See  Section  V.) 

•  Program  LSPRWC  will  issue  an  error  message  should  an  illegal  value  for  np  be 
entered.  (See  Section  V.) 

•  Program  LSPRWC  will  issue  an  error  message  should  an  illegal  value  for  nc  be 
entered.  (See  Section  V.) 

•  Program  LSPRWC  will  issue  an  error  message  should  the  number  of  (x,  y )  data 
pairs  be  insufficient  for  the  number  of  constraints,  £,  and  the  polynomial  order,  np. 
(See  Section  V.) 


8 


REFERENCES 


1.  Bovet,  D.  and  Cesati,  M.,  Understanding  the  LINUX  Kernel,  3rd  Edition,  ISBN-10: 
0-596-00565-2,  O’Reilly  Media,  Inc.,  2005. 

2.  Rosen,  K.  et  al.,  UNIX:  The  Complete  Reference,  2nd  Edition,  ISBN-10: 
0-07-226336-9,  MCGRAW-HILL,  2007. 

3.  Wikipedia,  “Least  squares,”  http://en.wikipedia.org/wiki/Least_squares/ 

4.  FORTRAN  77  Language  Reference  Manual.  Document  Number  007-0710-030,  Silicon 
Graphics,  Inc.,  1990. 

5.  Young.  D.A..  The  X  Window  System:  Programming  and  Applications  with  Xt, 
OSF/Motif.  2nd  Edition,  ISBN  0-13-123803-5,  Prentice-Hall,  1994. 

6.  http://gcc.gnu.org/fortran/ 

7.  Wikipedia,  “Mac  OS  X,”  http://en.wikipedia.org/wiki/Mac_OS_X 

8.  Wikipedia,  “Yellow  Dog  Linux,”  http://en.wikipedia.org/wiki/Yellow_Dog_Linux 

9.  Adams,  J.C.  et  al.,  FORTRAN  95  HANDBOOK.  ISBN  0-262-51096-0,  The  MIT  Press, 
1997. 

10.  “Intel  Fortran  Compiler  for  Linux,”  http://software.intel.com/en-us/articles/intel- 
fortrancompiler-for-linux-9x-manuals/ 

1 1 .  Weinstock,  R.,  Calculus  of  Variations  with  Applications  to  Physics  and  Engineering, 
New  York,  Dover  Publications,  Inc.,  pp.  6,  1974. 

12.  Carnahan,  B.;  Luther,  H.A.;  and  Wilkes,  J.L.,  Applied  Numerical  Methods,  New  York, 
John  Wiley  and  Sons,  pp.  281-96,  1969. 


9 


APPENDIX  A 
DERIVATION 


A  set  of  (x,y)  data  pairs  are  to  be  “curve  fit”  using  the  method  of  least  squares;  that  is,  an 
np  order  polynomial  is  chosen  such  that  the  sum  of  the  squared  errors  is  minimized  over  the  set 
of  (x,y)  data  pairs.  In  addition,  the  np  order  polynomial  must  satisfy  a  set  of  constraints; 
namely,  that  the  nc  polynomial  derivative  be  specified  at  a  given  x  location. 

Then  with  the  np  order  polynomial  given  by 


the  squared  error  to  be  minimized  is 


Hx)  =  2  bJxJ> 

7=0 


1=1 


(1) 


(2) 


m  np 

=2><-I>W)2’ 


(3) 


subject  to  the  constraints 


Gl{xl)  =  Pk‘{xl)  =  Cl 


(4) 


F 

=  1 


f- 


bjXj  =Cn 

AU-W  J 


(5) 


for  /  =  1,2,3,-  •  •,£ 


where 


nP  =  polynomial  order, 

C  =  number  of  constraints, 
m  =  number  of  (XA)  data  pairs, 
xi  =  Ith  constraint  location, 
ki  =  Ith  constraint  order  (0  <  ki  <  np ) 

Ci  =  value  of  the  hh  constraint. 


A-l 


Using  the  method  of  undetennined  Lagrange  multipliers,  the  relation  to  be  minimized  is 


£ 

S*=S+J^X,G,  (6) 

1=1 

by  setting 


as* 

dbj 


0  for  =  0, 1,*  •  *,np. 


From  equations  (3),  (5),  and  (6), 


(V) 


= 


i=  1 


X )bJX ft 

7=0 


7! 


1=  1 


j=kl 


(j~ki)\ 


bjA 


i-ki 


(8) 


Then  using  equation  (7), 


where 


as* 

a^ 


m  np  E, 

I>2*-')(y<  -  2  Vi)  +  E  =  0 


i=  1 


7=0  /=1 


F 


U  ~ 


J\  J~kl 
(j-k,y.  1 


0 


,  J>k, 
,  J  <  ki 


Rewriting  equation  (9)  gives 


E2> 

i=  1  7=0 


jx\  j  + 


w  Y.hFu  -  ^ 


1=  1 


i=  1 


(9) 


(10) 


(11) 


for  J  —  0, 1,  • 

The  minimization  equation  (11)  provides  np  +  1  equations  in  the  np  +  1  unknown 
coefficients  bj  and  c  unknown  Lagrange  multipliers  /l/.  Constraint  equation  (5)  provides  the 
remaining  £  equations  in  the  unknown  bj.  Since  equations  (5)  and  (11)  are  linear  and  algebraic, 
they  may  be  solved  simultaneously  using  a  matrix  inversion  technique. 
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APPENDIX  B 

FORTRAN  77  CODE 


c* *********************************************************************  * 
C*  * 

c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

C*  REF:  CARNAHAN, B.,  LUTHER,  H.A.,  AND  WILKES,  J.O.:  APPLIED  NUMERICAL  * 
C*  METHODS,  NEW  YORK,  JOHN  WILEY  &  SONS,  1969,  PP.  571-584.  * 

* 

WE INSTOCK,  R. :  CALCULUS  OF  VARIATIONS  WITH  APPLICATIONS  TO  * 

PHYSICS  AND  ENGINEERING,  NEW  YORK,  DOVER  PUBLICATIONS,  INC.,  * 
1974,  P.  6. 


LEAST  SQUARES  POLYNOMIAL  REGRESSION  WITH  CONSTRAINTS  (LSPRWC) 
WRITTEN  BY:  C.D.  MIKKELSEN 
2  MAY  1988 

AERODYNAMICS  TECHNOLOGY  BRANCH  (AMSMI-RD-SS-AT) 

SYSTEMS  SIMULATION  AND  DEVELOPMENT  DIRECTORATE 
US  ARMY  MISSILE  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING  CENTER 

US  ARMY  MISSILE  COMMAND 
REDSTONE  ARSENAL,  ALABAMA  35898-5252 


C* 

C* 

C* 

C*  1974,  P.  6.  * 

C*  * 

C*  REVISION  DATE:  7  JULY  2011  * 

C*  * 

c* *********************************************************************  * 

c 

c* *********************************************************************  * 
C*  * 

C*  PROGRAM  LSPRWC  REQUIRES  THE  FOLLOWING  SUBPROGRAMS:  * 

C*  * 

C*  CBUS06  CDS  DFCUS  FACTRL  FDFOS  GJEMPS  LSPRWC  * 

C*  PES  PODS  POWER  QPS  XDPEDS  XYDIS  YNOUS  * 

C*  * 

c* *********************************************************************  * 

c 

c* *********************************************************************  * 

c* 

C*  PARAMETERS: 

C* 

C*  MXNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  1  (MXNRM1  =  MXNR  -  1) 

C* 

C*  LOGICAL  UNIT  DEFINITIONS: 

C* 

C*  UNIT  FILE 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

c* *********************************************************************  * 

c 

c 

c 

c 

c 


1  INPUT  (X, Y)  DATA  PAIR  SCRATCH  FILE 

2  (XC,NC,CV)  CONSTRAINT  SCRATCH  FILE 

3  EVALUATED  DATA  SCRATCH  FILE 

4  TEMPORARY  USE 

5  STANDARD  INPUT 

6  STANDARD  OUTPUT 


PROGRAM  MAIN 

IMPLICIT  REAL*8(A-H,0-Z) 

CHARACTER  CCV*1 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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o  o  o  o 


INTEGER  DUM2 

00000061 

c 

00000062 

PARAMETER  (MXNRM1=29) 

00000063 

c 

00000064 

COMMON/B/B ( 0 : MXNRM1 ) 

00000065 

COMMON/NCONST/NCONST 

00000066 

COMMON / NEVAL / NEVAL 

00000067 

COMMON /NP/NP 

00000068 

COMMON/ SDEV/ SDEV 

00000069 

c 

00000070 

c 

00000071 

C***********************************************************************00000072 

c* 

*00000073 

c* 

OPEN  THE  SCRATCH  FILES 

*00000074 

c* 

*00000075 

C***********************************************************************00000076 

c 

00000077 

c 

00000078 

OPEN (UNIT=1,STATUS=' SCRATCH' , FORM= ' UNFORMATTED ' ,ERR=116) 

00000079 

OPEN (UNIT=2,STATUS=' SCRATCH' , FORM= ' UNFORMATTED ' ,ERR=117) 

00000080 

OPEN (UNIT=3,STATUS=' SCRATCH' , FORM= ' UNFORMATTED ' ,ERR=118) 

00000081 

c 

00000082 

c 

00000083 

C***********************************************************************00000084 

c* 

*00000085 

c* 

WRITE  THE  PROGRAM  DESCRIPTION 

*00000086 

c* 

*00000087 

q*  **********************************************************************  qqqqqq$% 

c 

00000089 

c 

00000090 

WRITE (UNIT=6 , FMT=201 ) 

00000091 

WRITE (UNIT=6 , FMT=202 ) 

00000092 

READ (UNI T= 5 , FMT=203 ) CCV 

00000093 

c 

00000094 

c 

00000095 

C***********************************************************************00000096 

c* 

*00000097 

c* 

INPUT  THE  X-Y  DATA  POINTS  TO  BE  FITTED 

*00000098 

c* 

*00000099 

q***********************************************************************qqqqoiqq 

c 

00000101 

c 

00000102 

101 

CALL  XYDIS(*119, *115) 

00000103 

c 

00000104 

c 

00000105 

q***********************************************************************qqqqqIQ6 

c* 

*00000107 

c* 

INPUT  THE  CONSTRAINTS 

*00000108 

c* 

*00000109 

c***********************************************************************ooooono 

c 

00000111 

c 

00000112 

102 

NCONST=0 

00000113 

103 

WRITE (UNIT=6 , FMT=208 ) 

00000114 

CALL  YNOUS (*104, *105, *103) 

00000115 

104 

CALL  CDS(*120, *115) 

00000116 

00000117 

00000118 


***********************************************************************00000119 

*  *00000120 
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C*  INPUT  THE  POLYNOMIAL  ORDER  *00000121 

C*  *00000122 

C***********************************************************************00000123 
C  00000124 

C  00000125 

105  CALL  PODS (*115)  00000126 

C  00000127 

C  00000128 

C***********************************************************************00000129 
C*  *00000130 

C*  PERFORM  THE  LEAST  SQUARES  POLYNOMIAL  REGRESSION  WITH  CONSTRAINTS  *00000131 

C*  *00000132 

C***********************************************************************00000133 
C  00000134 

C  00000135 

CALL  LSPRWC(*121)  00000136 

WRITE (UNIT=6,FMT=209)  00000137 

WRITE (UNIT=6 , FMT=210) (I,B(I) ,I=0,NP)  00000138 

WRITE (UNIT=6 , FMT=211 ) NP, SDEV  00000139 

WRITE (UNIT=6 , FMT=202 )  00000140 

READ (UNI T= 5 , FMT=203 ) CCV  00000141 

C  00000142 

C  00000143 

C***********************************************************************00000144 
C*  *00000145 

C*  INPUT  THE  X  DATA  POINTS  FOR  EVALUATION  *00000146 

C*  OF  THE  LEAST  SQUARES  POLYNOMIAL  *00000147 

C*  *00000148 

C***********************************************************************00000149 
C  00000150 

C  00000151 

106  WRITE (UNIT=6 , FMT=212 )  00000152 

LCV2=0  00000153 

CALL  YNOUS (*107, *112, *106)  00000154 

107  CALL  XDPEDS(LCV2, *122, *115)  00000155 

C  00000156 

C  00000157 

C***********************************************************************00000158 
C*  *00000159 

C*  EVALUATE  THE  LEAST  SQUARES  POLYNOMIAL  AT  THE  PRESCRIBED  X  VALUES  *00000160 

C*  *00000161 

C***********************************************************************00000162 
C  00000163 

C  00000164 

108  CALL  PES ( *123 )  00000165 

109  WRITE (UNIT=6 , FMT=205 )  00000166 

CALL  YNOUS (*110, *112, *109)  00000167 

110  REWIND  3  00000168 

WRITE (UNIT=6 , FMT=206 )  00000169 

DO  111  1=1 , NEVAL  00000170 

READ (UNIT=3 , ERR=124 ) X, Y  00000171 

IF (MOD (I, 20) .NE.O)  GO  TO  111  00000172 

WRITE (UNIT=6 , FMT=202 )  00000173 

READ (UNIT=5,FMT=203) CCV  00000174 

WRITE (UNIT=6 , FMT=206 )  00000175 

111  WRITE (UNIT=6 , FMT=207 ) I , X, Y  00000176 

WRITE (UNIT=6 , FMT=202 )  00000177 

READ (UNIT=5,FMT=203) CCV  00000178 

C  00000179 

C  00000180 


108  CALL  PES (*123) 

109  WRITE (UNIT=6 , FMT=205 ) 

CALL  YNOUS (*110, *112, *109) 

110  REWIND  3 

WRITE (UNIT=6 , FMT=206 ) 

DO  111  1=1, NEVAL 
READ (UNIT=3 , ERR=124 ) X, Y 
IF(MOD(I,20) .NE.O)  GO  TO  111 
WRITE (UNIT=6 , FMT=202 ) 

READ (UNIT=5 , FMT=203 ) CCV 
WRITE (UNIT=6 , FMT=206 ) 

111  WRITE (UNIT=6 , FMT=207 ) I , X, Y 
WRITE (UNIT=6 , FMT=202 ) 

READ (UNIT=5 , FMT=203 ) CCV 
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C***********************************************************************00000181 

C*  *00000182 

C*  PLOT  THE  RESULTS  *00000183 

C*  *00000184 

C***********************************************************************00000185 
C  00000186 

C  00000187 

112  WRITE (UNIT=6,FMT=213)  00000188 

CALL  YNOUS (*113,*115,*112)  00000189 

113  CALL  QPS(*125)  00000190 

GO  TO  115  00000191 

C  00000192 

C  00000193 

C***********************************************************************00000194 
C*  *00000195 

C*  SAVE  RESULTS  ON  FORMATTED  DISK  FILE  *00000196 

C*  *00000197 

C***********************************************************************00000198 
C  00000199 

C  00000200 

114  CALL  FDFOS(LCV2, *126)  00000201 

C  00000202 

C  00000203 

C***********************************************************************00000204 
C*  *00000205 

C*  RESTART  PROCEDURE  *00000206 

C*  *00000207 

C***********************************************************************00000208 
C  00000209 

C  00000210 

115  WRITE (UNIT=6,FMT=214)  00000211 

READ (UNIT=5 , FMT=204 ) LCV1  00000212 

GO  TO  (101,104,105,107,113,114) ,LCV1  00000213 

STOP  00000214 

C  00000215 

C  00000216 

C***********************************************************************00000217 
C*  *00000218 

C*  ERROR  MESSAGES  *00000219 

C*  *00000220 

C***********************************************************************00000221 
C  00000222 

C  00000223 

116  WRITE (UNIT=6 , FMT=215 )  00000224 

STOP  00000225 

117  WRITE (UNIT=6 , FMT=216 )  00000226 

STOP  00000227 

118  WRITE (UNIT=6 , FMT=217 )  00000228 

STOP  00000229 

119  WRITE (UNIT=6 , FMT=218 )  00000230 

STOP  00000231 

120  WRITE (UNIT=6 , FMT=219 )  00000232 

STOP  00000233 

121  WRITE (UNIT=6 , FMT=220)  00000234 

GO  TO  128  00000235 

122  WRITE (UNIT=6 , FMT=221 )  00000236 

STOP  00000237 

123  WRITE (UNIT=6 , FMT=222 )  00000238 

STOP  00000239 

124  WRITE (UNIT=6 , FMT=223 )  00000240 
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STOP 

00000241 

125 

WRITE (UNIT=6 , FMT=224 ) 

00000242 

STOP 

00000243 

126 

WRITE (UNIT=6 , FMT=225 ) 

00000244 

STOP 

00000245 

127 

WRITE (UNIT=6 , FMT=226 ) 

00000246 

STOP 

00000247 

128 

WRITE (UNIT=6 , FMT=202 ) 

00000248 

READ (UNI T= 5 , FMT=203 ) CCV 

00000249 

GO  TO  115 

00000250 

C 

00000251 

C 

00000252 

C***********************************************************************00000253 

c* 

*00000254 

c* 

FORMAT  STATEMENTS 

*00000255 

c* 

*00000256 

c 

c 

201 


202 

203 

204 

205 

206 

207 

208 

209 

210 
211 

212 

213 

214 


215 

216 
217 


00000258 

00000259 

FORMAT (/, T4 PROGRAM  LSPRWC  IS  AN  INTERACTIVE  FORTRAN  PROGRAM  TO  P00000260 
>ERFORM  A  LEAST SQUARES  POLYNOMIAL  REGRESSION  WITH  CONTRAINTS;  00000261 
>THAT  IS,  A  SET  OF  X-Y  DATA' , /,' POINTS  IS  CURVE  FIT  WITH  AN  NP  ORDE00000262 
>R  POLYNOMIAL  OF  THE  FORM' , //,T17 , 'P(X)=B0+B1*X+B2*X**2+B3*X**3+. . .00000263 
>.+BNP*X**NP' ,//, 'WITH  ANY  POLYNOMIAL  DERIVATIVES,  ZERO  THROUGH  NP, 00000264 

>  SPECIFIED  AT  GIVEN  X ',/,' LOCATIONS .  THE  PROCEDURE  FOLLOWED  IS  TH00000265 
>E  METHOD  OF  LEAST  SQUARES  USING' ,/, 'UNDETERMINED  LAGRANGE  MULTIPLI00000266 
>ERS. ' ,//,T4, 'AS  AN  INTERACTIVE  PROGRAM,  LSPRWC  IS  SELF-EXPLANATORY00000267 

>  AND  PROMPTS  FOR  THE ',/,' NECESSARY  INFORMATION.  THE  X-Y  DATA  TO  B00000268 
>E  FITTED  MAY  BE  ENTERED  BY  NAMELIST ',/,' OR  READ  FROM  A  FORMATTED  D00000269 
>ISC  FILE.  PROGRAM  LSPRWC  WILL  ALSO  EVALUATE  THE ',/,' RESULTANT  LEA00000270 
>ST  SQUARES  POLYNOMIAL  AT  PRESCRIBED  VALUES  OF  X  WHICH,  AGAIN, ',/,' 0000027 1 
>MAY  BE  ENTERED  BY  NAMELIST  OR  READ  FROM  A  FORMATTED  DISK  FILE.') 

FORMAT (/, T19, '-  ENTER/RETURN  TO  CONTINUE  -') 

FORMAT (Al) 

FORMAT (II) 

FORMAT (/,' SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION? 

>') 


00000272 
00000273 
00000274 
00000275 
( Y/N) 00000276 
00000277 


FORMAT (/,T10, 'NO. ' ,T26, 'X' ,T49, 'Y' )  00000278 

FORMAT (T8, 15, 1P2D23. 13)  00000279 

FORMAT(/, 'ARE  CONSTRAINTS  DESIRED?  (Y/N)')  00000280 

FORMAT (T7 , 'LEAST  SQUARES  POLYNOMIAL' ,//,T19, 'P(X)=B(0)+B(1) *X+B(2)00000281 
>*X**2+. . . . +B (NP) *X**NP ' ,//,T22, 'I',T36,'B(I)',/)  00000282 

FORMAT ( '  ' ,T22,I1,1PD25.13)  00000283 

FORMAT (/, T4 ,' THE  STANDARD  DEVIATION  FOR  THIS  POLYNOMIAL  OF  ORDER  00000284 

>  ',11,'  IS  ' , 1PD12 . 5 )  00000285 

FORMAT (/,' SHOULD  DATA  POINTS  BE  PRESCRIBED  FOR  EVALUATION  OF  THE  L00000286 

>EAST  SQUARES' ,/, 'POLYNOMIALS?  (Y/N)')  00000287 

FORMAT (/,' SHOULD  THE  RESULTS  OF  THIS  RUN  BE  QUICK-PLOTTED?  (Y/N) ') 00000288 

>  00000289 
FORMAT (/, 'ENTER: ',/, '1,  TO  RESTART  THE  PROGRAM' ,/,' 2 ,  TO  CHANGE  TH00000290 

>E  CONSTRAINTS ',/,' 3 ,  TO  CHANGE  THE  POLYNOMIAL  ORDER',/, '4,  TO  CHAN00000291 
>GE  THE  PRESCRIBED  X  VALUES  FOR  EVALUATION  OF  THE  POLYNOMIAL ',/,' 5 , 00000292 

>  TO  PLOT  THE  RESULTS ',/,' 6 ,  TO  SAVE  THE  RESULTS  ON  A  FORMATTED  DIS00000293 

>C  FILE',/, '7,  TO  STOP')  00000294 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  OPEN  ERROR  ON  UNIT  1,  STATUS  =  00000295 
>" SCRATCH" ' )  00000296 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  OPEN  ERROR  ON  UNIT  2,  STATUS  =  00000297 
>" SCRATCH" ' )  00000298 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  OPEN  ERROR  ON  UNIT  3,  STATUS  =  00000299 
>" SCRATCH" ' )  00000300 
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218  FORMAT ( 'SUBROUTINE  XYDIS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

219  FORMAT ( 'SUBROUTINE  CDS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

220  FORMAT ( 'SUBROUTINE  LSPRWC  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

221  FORMAT ( 'SUBROUTINE  XDPEDS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

222  FORMAT ( 'SUBROUTINE  PES  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

223  FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  UNFORMATTED  READ  ERROR  ON$UNIT 

>3') 

224  FORMAT ( 'SUBROUTINE  QPS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

225  FORMAT ( 'SUBROUTINE  FDFOS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

226  FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

END 


00000301 

00000302 

00000303 

00000304 

00000305 

00000306 

00000307 

00000308 

00000309 

00000310 

00000311 
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SUBROUTINE  CDS ( * , * ) 


C 

C 

C* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c 

c 

c 


********************************************************************** 

* 

CONSTRAINT  DEFINITION  SUBROUTINE  (CDS)  * 

* 

REVISION  DATE:  1  JULY  2011  * 

* 

********************************************************************** 


********************************************************************** 


SUBROUTINE  CDS  DEFINES  CONSTRAINTS  FOR  PROGRAM  LSPRWC.  * 

* 

INPUT/OUTPUT  VARIABLES:  * 

* 

CV  =  VALUE  OF  THE  NC-TH  DERIVATIVE  FOR  THE  CONSTRAINT  AT  XC  * 

NC  =  CONSTRAINT  ORDER  * 

NCONST  =  NUMBER  OF  CONSTRAINTS  * 

XC  =  X-VALUE  FOR  THE  CONSTRAINT  * 

* 

PARAMETERS :  * 

* 

MXNAC  =  MAXIMUM  NUMBER  OF  ALLOWABLE  CONSTRAINTS  * 

* 

********************************************************************** 


IMPLICIT  REAL*8 (A-H,0-Z) 

CH ARACTE R  CCV*l,DFmt*3,File*80,Fmt*80 
INTEGER  DUM2 


C 

PARAMETER  (R_NaN=-999 .E+00) 
PARAMETER  (MXNAC=20) 

C 

COMMON / C V/ C V ( MXNAC ) 

COMMON / NC / NC ( MXNAC ) 
COMMON/NCONST/NCONST 
COMMON / XC / XC ( MXNAC ) 

C 

DATA  DFmt/ ' (*) '/ 

C 

NAMELIST/PARM/CV, NC , XC 
C 

c 


c*  * 

C*  INPUT  THE  CONSTRAINTS  * 

C*  * 

c 

c 

DO  101  1=1, MXNAC 
CV ( I ) =R_NaN 
NC ( I ) =-l 

101  XC ( I ) =R_NaN 

WRITE (UNIT=6 , FMT=201 ) 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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READ (UNI T= 5 , FMT=202 ) LCV1 

00000061 

IF (LCV1 . EQ . 2 )  GO  TO  104 

00000062 

c 

00000063 

102 

WRITE (UNIT=6 , FMT=203 )MXNAC 

00000064 

WRITE (UNIT=6 , FMT=204 ) 

00000065 

CALL  NDRUS ( ' PARM' , 4 , *102 , *121 ) 

00000066 

READ (UNIT=4 , NML=PARM, ERR=102 ) 

00000067 

CALL  NDRUSE(*102,*121) 

00000068 

DO  103  NCONST=MXNAC , 1,-1 

00000069 

IF (XC (NCONST) . NE.R  NaN)  GO  TO  111 

00000070 

103 

CONTINUE 

00000071 

GO  TO  111 

00000072 

C 

00000073 

104 

WRITE (UNIT=6 , FMT=205 ) MXNAC 

00000074 

WRITE (UNIT=6 , FMT=206 ) 

00000075 

READ (UNIT=5 , FMT=207 ) FILE 

00000076 

105 

WRITE (UNIT=6 , FMT=208 ) 

00000077 

WRITE (UNIT=6 , FMT=209 ) 

00000078 

CALL  DFCUS ( Fmt , DFmt ,*105) 

00000079 

OPEN (UNIT=4 , FILE=File, STATUS= ' OLD ' , ERR=122 ) 

00000080 

REWIND  4 

00000081 

IF ( Fmt (1:3) . EQ . DFmt )  GO  TO  107 

00000082 

DO  106  NCONST=l, MXNAC 

00000083 

106 

READ (UNIT=4 , FMT=Fmt , END=109 , ERR=123 ) XC (NCONST) 

, NC (NCONST) ,  00000084 

>CV( NCONST) 

00000085 

GO  TO  110 

00000086 

107 

DO  108  NCONST=l, MXNAC 

00000087 

108 

READ (UNIT=4,FMT=*,END=109,ERR=123)XC (NCONST) , NC (NCONST) ,CV (NCONST) 00000088 

GO  TO  110 

00000089 

109 

NCONST=NCONST-l 

00000090 

110 

CLOSE (UNIT=4 , STATUS= ' KEEP ' , ERR=124 ) 

00000091 

C 

00000092 

C 

00000093 

q*  **********************************************************************  qqqqqq^  ^ 

c* 

*00000095 

c* 

SORT  AND  EDIT  THE  CONSTRAINTS 

*00000096 

c* 

*00000097 

Q*  ********************************************************************** QQQQQQ$% 

c 

00000099 

c 

00000100 

111 

IF ( NCONST. EQ.O)  GO  TO  120 

00000101 

IF ( NCONST. EQ.l)  GO  TO  117 

00000102 

DO  112  1=1 , NCONST-1 

00000103 

DO  112  J=I+1, NCONST 

00000104 

IF(XC(I) .NE.XC(J) )  GO  TO  112 

00000105 

IF(NC(I) .NE.NC(J) )  GO  TO  112 

00000106 

XC ( I ) =R  NaN 

00000107 

112 

CONTINUE 

00000108 

113 

LCV1=0 

00000109 

DO  115  1=2, NCONST 

00000110 

IF (XC ( 1-1 ) . LT . XC ( I ) )  GO  TO  115 

00000111 

IF (XC ( 1-1 ) . GT . XC ( I ) )  GO  TO  114 

00000112 

IF  (NC  ( 1-1 )  .  LT .  NC  ( I ) )  GO  TO  115 

00000113 

114 

LCV1=1 

00000114 

DUM1=XC ( 1-1 ) 

00000115 

DUM2=NC ( 1-1 ) 

00000116 

DUM3=CV ( 1-1 ) 

00000117 

XC(I-1)=XC(I) 

00000118 

NC(I-1)=NC(I) 

00000119 

CV(I-1)=CV(I) 

00000120 

B-8 


XC ( I ) =DUM1 

00000121 

NC ( I ) =DUM2 

00000122 

CV ( I ) =DUM3 

00000123 

115 

CONTINUE 

00000124 

IF(LCVl.EQ.l)  GO  TO  113 

00000125 

DO  116  NCONST=MXNAC ,1,-1 

00000126 

IF (XC (NCONST) . NE.R  NaN)  GO  TO  117 

00000127 

116 

CONTINUE 

00000128 

C 

00000129 

117 

WRITE (UNIT=6 , FMT=210) 

00000130 

CALL  YNOUS (*118, *120, *117) 

00000131 

118 

WRITE (UNIT=6 , FMT=211 ) 

00000132 

DO  119  1=1, NCONST 

00000133 

IF (MOD (I, 20) .NE.O)  GO  TO  119 

00000134 

WRITE (UNIT=6 , FMT=212 ) 

00000135 

READ (UNIT=5 , FMT=213 ) CCV 

00000136 

WRITE ( UNIT=6 , FMT=2 1 1 ) 

00000137 

119 

WRITE (UNIT=6 , FMT=214 ) I , XC ( I ) ,NC(I) ,CV(I) 

00000138 

WRITE (UNIT=6 , FMT=212 ) 

00000139 

RE AD ( UN I T= 5 , FMT= 2 1 3 ) CCV 

00000140 

120 

RETURN 

00000141 

C 

00000142 

C 

00000143 

C***********************************************************************00000144 

c* 

*00000145 

c* 

ERROR  MESSAGES 

*00000146 

c* 

*00000147 

C***********************************************************************00000148 

c 

00000149 

c 

00000150 

121 

WRITE (UNIT=6 , FMT=215 ) 

00000151 

GO  TO  126 

00000152 

122 

WRITE (UNIT=6 , FMT=216 ) 

00000153 

CALL  CBUS06 (6 , File, *125) 

00000154 

RETURN  1 

00000155 

123 

WRITE (UNIT=6 , FMT=217 ) 

00000156 

CALL  CBUS06 (6 , File, *125) 

00000157 

RETURN  1 

00000158 

124 

WRITE (UNIT=6 , FMT=218 ) 

00000159 

CALL  CBUS06 (6 , File, *125) 

00000160 

RETURN  1 

00000161 

125 

WRITE (UNIT=6 , FMT=219 ) 

00000162 

RETURN  1 

00000163 

126 

WRITE (UNIT=6 , FMT=212 ) 

00000164 

RE AD ( UN I T= 5 , FMT= 2 1 3 ) CCV 

00000165 

RETURN  2 

00000166 

c 

c 


00000167 

00000168 


C***********************************************************************00000169 

C*  *00000170 

C*  FORMAT  STATEMENTS  *00000171 

C*  *00000172 

C***********************************************************************00000173 

C  00000174 

C  00000175 

201  FORMAT (/,' SELECT  THE  SOURCE  OF  INPUT  FOR  THE  CONSTRAINTS  FROM  THE  00000176 

FOLLOWING  LIST: ' ,11, '1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 2 ,  FOR00000177 
>  FORMATTED  DISK  FILE  INPUT ',///' NOTE :  THE  CONSTRAINTS  ARE  COMPLETE00000178 

>LY  INDEPENDENT  OF  THE  X-Y  DATA  POINTS  TO  BE FITTED .' )  00000179 

202  FORMAT ( II )  00000180 


B-9 


203 


204 

205 


206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 

218 
219 


FORMAT (/,' INPUT  THE  CONSTRAINTS  (',12,'  MAX)  BY  NAMELIST  WHERE :', /00000181 
>/,'CV  =  ARRAY  OF  POLYNOMIAL  DERIVATIVE  VALUES ',/,' NC  =  ARR00000182 
>AY  OF  POLYNOMIAL  DERIVATIVE  ORDERS  (CONSTRAINT  ORDER) ',/,' XC  =00000183 
>  ARRAY  OF  POLYNOMIAL  DERIVATIVE  LOCATIONS  (CONSTRAINT  X  VALUES) ', /00000184 
>/ , ' NOTE :  NC  VALUES  MUST  BE  IN  THE  RANGE  0  TO  NP .',//, T2 ,' CURRENT  V00000185 
>ALUES  ARE : ' )  00000186 


FORMAT  (  /  ,  '  $PARM  CV= _ , _ , _ ,NC= _ , _ , _ ,XC= _ , _ ,_00000187 

> _ $',/)  00000188 

FORMAT (/,' CONSTRAINTS  (',12,'  MAX)  ARE  INPUT  FROM  FORMATTED  DISK  F00000189 
>ILES  AS  (XC,NC,CV)  TRIPLES ',/, 'WHERE CV  =  POLYNOMIAL  DER00000190 
>IVATIVE  VALUE',/, 'NC  =  POLYNOMIAL  DERIVATIVE  OR  CONSTRAINT  0RD00000191 
>ER',/,'XC  =  POLYNOMIAL  DERIVATIVE  LOCATION  (CONSTRAINT  X  VALUE00000192 

>)',//, 'NOTE:  NC  VALUES  MUST  BE  IN  THE  RANGE  0  TO  NP.')  00000193 

FORMAT (/,' INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE  DATA  SET. 00000194 

00000195 
00000196 
00000197 
F000000198 
00000199 
(Y00000200 
00000201 
00000202 
00000203 
00000204 
00000205 
00000206 
=  00000207 
00000208 


>') 

FORMAT (A80) 

FORMAT (/,' INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES)') 

FORMAT (' NOTE :  AN  EXAMPLE  FORMAT  IS  " (2E15 . 6 )"',/, T7 ,' ENTER" (*) " 

>R  A  FREE  FIELD  READ  (DEFAULT  FORMAT)') 

FORMAT (/,' SHOULD  THE  CONSTRAINTS  BE  DISPLAYED  FOR  VERIFICATION? 

>/N) ' ) 

FORMAT (/,T6, 'NO. ' ,T27, ' XC ( I ) ' , T35 , ' NC ( I ) ' ,T58, 'CV(I) ' ) 

FORMAT (/, T19, '-  ENTER/RETURN  TO  CONTINUE  -') 

FORMAT (Al) 

FORMAT (T8 , 13 , 1PD23 .13, 17, D2 3. 13) 

FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 

FORMAT ( / , ' ERROR  IN  PROGRAM  LSPRWC :  OPEN  ERROR  ON  UNIT  4 ,  STATUS 
>"OLD" ,  FILE  =' ,/) 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  READ  ERROR  ON  UNIT  4,  FILE  =',/00000209 
>)  00000210 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  CLOSE  ERROR  ON  UNIT  4,  FILE  =',00000211 

>/)  00000212 

FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  PROGRAM  LSPRWC')  00000213 

END  00000214 
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FUNCTION  FACTRL(N) 

00000001 

c 

00000002 

c 

00000003 

C***********************************************************************00000004 

c* 

*00000005 

c* 

FACTORIAL  EVALUATION  FUNCTION  SUBPROGRAM  (FACTRL) 

*00000006 

c* 

*00000007 

c* 

REVISION  DATE:  15  JUNE  1989 

*00000008 

c* 

*00000009 

c***********************************************************************ooooooio 

c 

00000011 

c 

00000012 

INTEGER  FACTRL 

00000013 

FACTRL=1 

00000014 

IF(N.LE.l)  RETURN 

00000015 

DO  101  1=0, N-l 

00000016 

101 

FACTRL=FACTRL* (N-I) 

00000017 

END 

00000018 

B-ll 


SUBROUTINE  FDFOS (LCV2 , * )  00000001 
C  00000002 
C  00000003 
C***********************************************************************00000004 
C*  *00000005 
C*  FORMATTED  DISK  FILE  OUTPUT  SUBROUTINE  (FDFOS)  *00000006 
C*  *00000007 
C*  REVISION  DATE:  2  OCTOBER  2007  *00000008 
C*  *00000009 
C***********************************************************************00000010 
C  00000011 
C***********************************************************************00000012 
C*  *00000013 
C*  PARAMETERS:  *00000014 
C*  *00000015 
C*  MXANC  =  MAXIMUM  NUMBER  OF  ALLOWABLE  CONSTRAINTS  *00000016 
C*  MXNR  =  MAXIMUM  NUMBER  OF  ROWS  *00000017 
C*  MXNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  1  (MXNRM1  =  MXNR  -  1)  *00000018 
C*  *00000019 
C***********************************************************************00000020 


c 

c 

c 

c 

c 


IMPLICIT  REAL*8(A-H,0-Z) 

CHARACTER  Alist*l,Blank*l,DFmt*9,File*80,Fmt*80,SList*l 

PARAMETER  (MXANC=20,MXNRM1=29) 

COMMON/B/B ( 0 : MXNRM1 ) 

COMMON /CV/CV ( MXANC ) 

COMMON / NC / NC ( MXANC ) 

COMMON/NCONST/NCONST 
COMMON / ND AT A/ ND AT A 
COMMON / NEVAL / NEVAL 
COMMON /NP/NP 
COMMON/ SDEV/ SDEV 
COMMON / XC / XC ( MXANC ) 

DIMENSION  AList(13) ,NList(13) , SList ( -1 : 11 ) , YK( -1 : 11 ) 
SAVE  Blank, DFmt 


00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 


00000042 

DATA  Blank/'  ' / ,DFmt/ ' ( 13E15 . 6) ' / , SList/"X" , "0" , " 1" , "2" , "3 " ,  "4" ,  00000043 

>"5" , "6" , "7" , "8" , "9" , "R" , "S"/  00000044 

C  00000045 

NAMELIST/PARM/AList  00000046 

C  00000047 

C  00000048 

C***********************************************************************00000049 
C*  *00000050 

C*  SELECT  AN  OUTPUT  OPTION  *00000051 

C*  *00000052 

C***********************************************************************00000053 
C  00000054 

C  00000055 

WRITE (UNIT=6 , FMT=201 )  00000056 

READ (UNIT=5 , FMT=202 ) FILE  00000057 

WRITE (UNIT=6 , FMT=203 )  00000058 

READ (UNIT=5 , FMT=204 ) LCV1  00000059 

IF (LCV1 . EQ . 2 .AND . LCV2 . EQ. 0)  GO  TO  113  00000060 
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c 

c 


OPEN (UNIT=4,FILE=File,F0RM=' FORMATTED' , STATUS= ' UNKNOWN ' ,ERR=114)  00000061 

REWIND  4  00000062 

IF (LCV1 . EQ . 2 )  GO  TO  104  00000063 

00000064 

00000065 


C***********************************************************************00000066 

C*  *00000067 

C*  OUTPUT  ALL  INPUT  DATA  AND  RESULTS  *00000068 

C*  *00000069 

C***********************************************************************00000070 


c 

c 

WRITE (UNIT=4 , FMT=205 , ERR= 
WRITE (UNIT=4 , FMT=206 , ERR= 
REWIND  1 

DO  101  1=1 , NDATA 
READ (UNIT=1 , ERR=116 ) X, Y 

101  WRITE (UNIT=4 , FMT=207 , ERR= 
IF ( NCONST . EQ . 0 )  GO  TO  102 
WRITE (UNIT=4 , FMT=208 , ERR= 
WRITE (UNIT=4 , FMT=209 , ERR= 

102  WRITE (UNIT=4 , FMT=210 , ERR= 
WRITE ( UNIT=4 , FMT=2 1 1 , ERR= 
WRITE (UNIT=4 , FMT=212 , ERR= 
IF (LCV2 . EQ . 0 )  GO  TO  112 
WRITE (UNIT=4 , FMT=213 , ERR= 
REWIND  3 

DO  103  1=1 , NEVAL 
READ (UNIT=3 , ERR=117 ) X,  Y 

103  WRITE (UNIT=4 , FMT=207 , ERR= 
GO  TO  112 

C 

C 

C*  ****************************  * 


115) 

00000071 

00000072 

00000073 

115) 

00000074 

115)1, X,Y 

00000075 

00000076 

00000077 

00000078 

115) 

00000079 

00000080 

115) (I,XC(I) ,NC(I) ,CV(I) ,1=1, NCONST) 

00000081 

115) 

00000082 

115) (I,B(I) , 1=0 ,NP ) 

00000083 

115)NP,SDEV 

00000084 

115) 

00000085 

00000086 

115)1, X,Y 

00000087 

00000088 

00000089 

00000090 

00000091 

00000092 

00000093 

*****************************************00000094 

c* 

C*  OUTPUT  ONLY  THE  EVALUATED  DATA 

C* 

C*  **************************************************  * 

c 

c 

104  DO  105  K=3 , 13 
AList (K) =Blank 

105  NList (K) =-999 
AList (1)='X' 

AList (2 ) = ' 0 ' 

WRITE (UNIT=6 , FMT=214 ) SList ( -1 ) 

WRITE (UNIT=6 , FMT=215 ) (SList(K) ,K,K=0,NP) 

WRITE (UNIT=6 , FMT=216 ) (SList(K) ,K=10, 11) 

WRITE (UNIT=6 , FMT=217 ) 

CALL  NDRUS ( ' PARM' ,9, *104, *118) 

READ (UNIT=9 , NML=PARM, ERR=104 ) 

CALL  NDRUSE(*104, *118) 

DO  106  NL=1 , 13 

IF(AList(NL) .EQ. Blank)  GO  TO  107 

106  CONTINUE 

107  NL=NL-1 
IF(NL.EQ.O)  GO  TO  119 
DO  110  1=1, NL 

DO  108  K=-l , 11 

IF(AList(I) .EQ.SList(K) )  GO  TO  109 


*00000095 

*00000096 

*00000097 

*******************00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

00000119 

00000120 
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108 

CONTINUE 

00000121 

GO  TO  120 

00000122 

109 

NList ( I ) =K 

00000123 

IF(K.EQ.IO)  NList ( I ) =NP+1 

00000124 

IF(K.EQ.ll)  NList ( I ) =NP+2 

00000125 

IF(NP.LT.K.AND.K.LT.IO)  GO  TO  120 

00000126 

110 

CONTINUE 

00000127 

C 

00000128 

WRITE (UNIT=6 , FMT=218 ) 

00000129 

CALL  CBUS06 ( 6 , DFmt , * 122 ) 

00000130 

CALL  DFCUS ( Fmt , DFmt ,*104) 

00000131 

REWIND  3 

00000132 

DO  111  1=1 , NEVAL 

00000133 

READ (UNIT=3 , ERR=117 ) ( YK ( K ) , K=- 1 , NP+2 ) 

00000134 

111 

WRITE (UNIT=4 , FMT=Fmt , ERR=115 ) ( YK(NList (K) ) ,K=1,NL) 

00000135 

112 

CLOSE (UNIT=4 , STATUS= ' KEEP ' , ERR=121 ) 

00000136 

RETURN 

00000137 

C 

00000138 

C 

00000139 

C***********************************************************************00000140 

c* 

*00000141 

c* 

ERROR  MESSAGES 

*00000142 

c* 

*00000143 

C***********************************************************************00000144 

c 

00000145 

c 

00000146 

113 

WRITE (UNIT=6 , FMT=219 ) 

00000147 

CALL  CBUS06 (6 , File, *122) 

00000148 

WRITE (UNIT=6 , FMT=220 ) 

00000149 

RETURN 

00000150 

114 

WRITE (UNIT=6 , FMT=221 ) 

00000151 

CALL  CBUS06 (6 , File, *122) 

00000152 

RETURN  1 

00000153 

115 

WRITE (UNIT=6 , FMT=222 ) 

00000154 

CALL  CBUS06 (6, File, *122) 

00000155 

RETURN  1 

00000156 

116 

WRITE (UNIT=6 , FMT=223 ) 

00000157 

RETURN  1 

00000158 

117 

WRITE (UNIT=6 , FMT=224 ) 

00000159 

RETURN  1 

00000160 

118 

WRITE (UNIT=6 , FMT=225 ) 

00000161 

RETURN  1 

00000162 

119 

WRITE (UNIT=6 , FMT=226 ) 

00000163 

GO  TO  104 

00000164 

120 

WRITE (UNIT=6 , FMT=227 ) I , AList ( I ) 

00000165 

GO  TO  104 

00000166 

121 

WRITE (UNIT=6 , FMT=228 ) 

00000167 

CALL  CBUS06 (6, File, *122) 

00000168 

RETURN  1 

00000169 

122 

WRITE (UNIT=6 , FMT=229 ) 

00000170 

RETURN  1 

00000171 

C 

00000172 

C 

00000173 

C***********************************************************************00000174 

c* 

*00000175 

c* 

FORMAT  STATEMENTS 

*00000176 

c* 

*00000177 

C***********************************************************************00000178 

C  00000179 

C  00000180 
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201 

202 

203 


204 

205 


206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 


217 

218 


219 


220 

221 

222 

223 

224 

225 

226 

227 

228 
229 


FORMAT (/,' INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE.')  00000181 

FORMAT (A80)  00000182 

FORMAT (/,' SELECT  THE  TYPE  OF  OUTPUT  FOR  THE  FORMATTED  DISK  FILE  FR00000183 
>OM  THE  FOLLOWING  LIST:',/,'l,  FOR  ALL  INPUT  DATA  AND  RESULTS 200000184 
>,  FOR  ONLY  THE  EVALUATED  DATA')  00000185 

FORMAT ( II )  00000186 

FORMAT (T6, 'LEAST  SQUARES  POLYNOMIAL  REGRESSION  WITH  CONSTRAINTS  (L00000187 
>SPRWC) ' ,//,T22, 'WRITTEN  BY:  C.D.  MIKKELSEN' , //,T31, '2  MAY  1988 ' , //00000188 
>,T13, 'AERODYNAMICS  TECHNOLOGY  BRANCH  (AMSMI-RD-SS-AT) ' , / , T13 , ' SYST00000189 
>EMS  SIMULATION  AND  DEVELOPMENT  DIRECTORATE ',/, T6 , 'US  ARMY  MISSILE  00000190 
> RESEARCH,  ENGINEERING,  AND  DEVELOPMENT  CENTER' , /,T25, 'US  ARMY  MISS00000191 
>ILE  COMMAND' ,/,T 18, 'REDSTONE  ARSENAL,  ALABAMA  35898-5252')  00000192 

FORMAT (/, T20, 'X-Y  DATA  POINTS  TO  BE  FITTED :',//, T12 ,' NO .', T28 ,' X ', 00000193 
>T51 , ' Y ' )  00000194 

FORMAT (T10, 15, 1P2D23. 13)  00000195 

FORMAT ( / , T25 , 'POLYNOMIAL  CONSTRAINTS' ,//,T8, 'NO. ' ,T29, 'XC(I) ' , T37 , 00000196 
>'NC(I) ' ,T60, 'CV(I) ' )  00000197 

FORMAT (T8 , 13 , 1PD23 . 13 , 17 ,D23 . 13 )  00000198 

FORMAT ( / , T28 , 'LEAST  SQUARES  POLYNOMIAL' ,//,T15, 'P(X)=B(0)+B(1) *X+B00000199 
>(2) *X**2+. . . . +B(NP) *x**NP' , // , T23 , 'I' ,T37, 'B(I) ' ,/)  00000200 

FORMAT (T2 3 , II , 1PD25 . 13 )  00000201 

FORMAT (/,T1, 'THE  STANDARD  DEVIATION  FOR  THIS  POLYNOMIAL  OF  ORDER  '00000202 
>,I1,'  IS  ' , 1PD12 . 5 )  00000203 

FORMAT (/,T22, 'X-Y  EVALUATED  DATA  POINTS :'//, T12 ,' NO .', T28 ,' X ', T51 , 00000204 
> ' Y ' )  00000205 

FORMAT (/,' INPUT  THE  EVALUATED  DATA  VARIABLE  LIST  FOR  OUTPUT  BY  NAM00000206 
>ELIST  WHERE: ,A1, ' "  =  X  VALUES')  00000207 

FORMAT ( ' " ' , A1 , ' "  =  ',11, 'th  DERIVATIVE  OF  Y  AT  X')  00000208 

FORMAT ( ' " ' , A1 , ' "  =  RADIUS  OF  CURVATURE  AT  X' ,A1, ' "  =  A00000209 

>RC  LENGTH' ,//, 'NOTE:  VARIABLES  WILL  BE  LISTED  IN  THE  ORDER  OF  INPU00000210 
>T ',//,' CURRENT  VALUES  ARE :')  00000211 

FORMAT (/, '$PARM  AList=  "X","0"$',/)  00000212 

FORMAT (/,' INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES) ',//, 'N000000213 
>TE :  AN  EXAMPLE  FORMAT  IS  " (F15 . 7 , E15 . 6 ) " ' , / , T7 , ' THE  DEFAULT  FORMAT00000214 
>  IS: ' )  00000215 

FORMAT (/, 'WARNING  IN  SUBROUTINE  FDFOS :  NO  DATA  POINTS  HAVE  BEEN  PR00000216 
>ESCRIBED  FOR  EVALUATION' ,/, 'OF  THE  LEAST  SQUARES  POLYNOMIAL.  HENCE00000217 
>,  FILE',/)  00000218 

FORMAT (/, 'WILL  NOT  BE  OPENED.')  00000219 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  OPEN  ERROR  ON  UNIT  4,  STATUS  00000220 
>=  "UNKNOWN",  FILE  =',/)  00000221 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  WRITE  ERROR  ON  UNIT  4,  FILE  =00000222 
>',/)  00000223 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  UNFORMATTED  READ  ERROR  ON  UNI00000224 
>T  1')  00000225 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  UNFORMATTED  READ  ERROR  ON  UNI00000226 
>T  3')  00000227 

FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  SUBROUTINE  FDFOS')  00000228 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  THE  OUTPUT  LIST  MUST  NOT  BE  E00000229 
>MPTY ' )  00000230 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  ILLEGAL  VALUE  FOR' ,/, 'AList (' 00000231 

>,  12 ,  ' )  =  ' , Al)  00000232 

FORMAT (/,' ERROR  IN  SUBROUTINE  FDFOS:  CLOSE  ERROR  ON  UNIT  4,  FILE  =00000233 

>',/)  00000234 

FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  SUBROUTINE  FDFOS')  00000235 

END  00000236 
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SUBROUTINE  GJEMPS ( EPS , INDIC , N , DETER, * ) 

C 

C 

C*  ********************************************************************  * 

c* 

C*  GAUSS-JORDAN  ELIMINATION  USING  MAXIMUM  PIVOT  STRATEGY  (GJEMPS) 

C* 

C*  REVISION  DATE:  21  APRIL  2000 

C* 

C*  ********************************************************************  * 

c 

c* ********************************************************************  * 

c* 

C*  SUBROUTINE  GJEMPS  COMPUTES  THE  INVERSE  OF  AN  N  BY  N  MATRIX  AND/OR 
C*  THE  N  SOLUTIONS  OF  A  SET  OF  N  LINEAR  EQUATIONS  IN  N  UNKNOWNS  USING 
C*  GAUSS-JORDAN  ELIMINATION  WITH  MAXIMUM  PIVOT  STRATEGY. 

C* 

C*  INPUT  VARIABLES: 

C* 

C*  EPS  =  MINIMUM  ALLOWABLE  MAGNITUDE  FOR  A  PIVOT  ELEMENT  (SUGGESTED 
C*  VALUE  1.0E-10) 

C*  INDIC  =  CONTROL  VARIABLE  SUCH  THAT: 

=  NEGATIVE,  TO  COMPUTE  THE  INVERSE  OF  N  BY  N  MATRIX  A  IN 
PLACE 

=  ZERO  ,  TO  COMPUTE  THE  N  SOLUTIONS  X(1)....X(N) 
CORRESPONDING  TO  THE  SET  OF  LINEAR  EQUATIONS  WITH 
AUGMENTED  MATRIX  OF  COEFFICIENTS  IN  THE  N  BY  N+l  ARRAY  A 
AND  IN  ADDITION  COMPUTE  THE  INVERSE  OF  THE  COEFFICIENT 
MATRIX  IN  PLACE 

=  POSITIVE,  TO  COMPUTE  THE  N  SOLUTIONS  X(1)....X(N) 
CORRESPONDING  TO  THE  SET  OF  LINEAR  EQUATIONS  WITH 
AUGMENTED  MATRIX  OF  COEFFICIENTS  IN  THE  N  BY  N+l  ARRAY  A 
=  NUMBER  OF  ROWS  IN  MATRIX  A 


C* 

C* 

c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

C*  N 
C* 

C*  OUTPUT  VARIABLES: 

C* 

C*  DETER  =  DETERMINATE  OF  THE  ORIGINAL  COEFFICIENT  MATRIX  FORMED  BY 
C*  THE  FIRST  N  COLUMNS  OF  ARRAY  A 

C*  X  =  VECTOR  OF  N  SOLUTIONS 

C* 

C*  INPUT/OUTPUT  VARIABLES: 

C* 

C*  A  =  AUGMENTED  COEFFICIENT  MATRIX 

C* 

C*  PARAMETERS: 

C* 

C*  MXNC  =  MAXIMUM  NUMBER  OF  COLUMNS  (MXNC  =  MXNR  +  1) 

C*  MXNCM1  =  MAXIMUM  NUMBER  OF  COLUMNS  MINUS  1  (MXNCM1  =  MXNC  -  1) 

C*  MXNR  =  MAXIMUM  NUMBER  OF  ROWS 

C*  MXNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  1  (MXNRM1  =  MXNR  -  1) 

C* 

C*  NOTE:  SHOULD  NO  ACCEPTABLE  PIVOT  ELEMENT  BE  FOUND,  COMPUTATIONS  ARE 
C*  TERMINATED  AND  THE  DETERMINATE  IS  RETURNED  WITH  A  TRUE  ZERO 

C*  VALUE. 

C* 

C*  REF:  CARNAHAN,  B.,  LUTHER,  H.A.,  AND  WILKES,  J.O.:  APPLIED 
C*  NUMERICAL  METHODS,  NEW  YORK,  1969,  JOHN  WILEY  &  SONS,  INC. 

C* 

C*  ********************************************************************  * 

c 


00000001 

00000002 

00000003 

*00000004 

*00000005 

*00000006 

*00000007 

*00000008 

*00000009 

*00000010 

00000011 

*00000012 

*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 

*00000025 

*00000026 

*00000027 

*00000028 

*00000029 

*00000030 

*00000031 

*00000032 

*00000033 

*00000034 

*00000035 

*00000036 

*00000037 

*00000038 

*00000039 

*00000040 

*00000041 

*00000042 

*00000043 

*00000044 

*00000045 

*00000046 

*00000047 

*00000048 

*00000049 

*00000050 

*00000051 

*00000052 

*00000053 

*00000054 

*00000055 

*00000056 

*00000057 

*00000058 

*00000059 

00000060 
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o  o 


C  00000061 

IMPLICIT  REAL*8(A-H,0-Z)  00000062 

C  00000063 

PARAMETER  (MXNC=31,MXNCM1=30,MXNR=30,MXNRM1=29)  00000064 

C  00000065 

DIMENSION  A(MXNR,MXNC) , IROW(MXNR) , JCOL (MXNR) , JORD (MXNR) ,X(MXNR) ,  00000066 

>Y (MXNR)  00000067 

C  00000068 

COMMON /B /B ( 0 : MXNRM1 )  00000069 

COMMON/C/C ( 0 : MXNRM1 , 0 : MXNCM1 )  0000007  0 

C  00000071 

EQUIVALENCE  (A(l, 1) ,C(0,0) ) , (X(l) ,B(0) )  00000072 

C  00000073 

C  00000074 

C***********************************************************************00000075 
C*  *00000076 

C*  BEGIN  ELIMINATION  PROCEDURE  *00000077 

C*  *00000078 

C***********************************************************************00000079 


c 

00000080 

c 

00000081 

IF(N.GT.MXNR)  GO  TO  114 

00000082 

MAX=N 

00000083 

IF ( INDIC . GE . 0 )  MAX=N+1 

00000084 

DETER=1 . 0D+00 

00000085 

DO  107  K=1 , N 

00000086 

KM1=K-1 

00000087 

c 

00000088 

c 

00000089 

C***********************************************************************00000090 

C*  *00000091 

C*  SEARCH  FOR  THE  PIVOT  ELEMENT  *00000092 

C*  *00000093 

C***********************************************************************00000094 
C  00000095 

C  00000096 

PIVOT=0 . 0D+00  00000097 

DO  103  1=1, N  00000098 

DO  103  J=1,N  00000099 

C  00000100 

C  00000101 

C***********************************************************************00000102 
C*  *00000103 

C*  SCAN  IROW  AND  JCOL  ARRAYS  FOR  INVALID  PIVOT  SUBSCRIPTS  *00000104 

C*  *00000105 

C***********************************************************************00000106 
C  00000107 

C  00000108 

IF(K.EQ.l)  GO  TO  102  00000109 

DO  101  ISCAN=1 , KM1  00000110 

DO  101  JSCAN=1 , KM1  00000111 

IF ( I . EQ . IROW ( ISCAN) )  GO  TO  103  00000112 

101  IF ( J. EQ. JCOL ( JSCAN) )  GO  TO  103  00000113 

102  IF (DABS ( A( I , J) ) .LE. DABS (PIVOT) )  GO  TO  103  00000114 

PIVOT=A( I , J)  00000115 

IROW(K) =1  00000116 

JCOL (K) =J  00000117 

103  CONTINUE  00000118 

00000119 
00000120 
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C***********************************************************************00000121 

C*  *00000122 

C*  INSURE  THE  SELECTED  PIVOT  IS  LARGER  THAN  EPS  *00000123 

C*  *00000124 

C***********************************************************************00000125 
C  00000126 

C  00000127 

IF ( DABS ( PIVOT ) .GT. EPS)  GO  TO  104  00000128 

DETER=0 . 0D+00  00000129 

RETURN  00000130 

C  00000131 

C  00000132 

C***********************************************************************00000133 
C*  *00000134 

C*  UPDATE  THE  DETERMINATE  VALUE  AND  NORMALIZE  PIVOT  ROW  ELEMENTS  *00000135 
C*  *00000136 

C***********************************************************************00000137 
C  00000138 

C  00000139 

104  DETER=DETER* PIVOT  00000140 

DO  105  J=1 , MAX  00000141 

105  A( IROW(K) , J) =A( IROW(K) , J) /PIVOT  00000142 

C  00000143 

C  00000144 

C***********************************************************************00000145 
C*  *00000146 

C*  CARRY  OUT  ELIMINATION  AND  DEVELOP  INVERSE  *00000147 

C*  *00000148 

C***********************************************************************00000149 
C  00000150 

C  00000151 

A( IROW(K) , JCOL(K) )=1 . OD+OO/PIVOT  00000152 

DO  107  1=1, N  00000153 

AI JCK=A( I , JCOL (K) )  00000154 

IF ( I . EQ. IROW(K) )  GO  TO  107  00000155 

A( I , JCOL (K) ) =-AI JCK/PIVOT  00000156 

DO  106  J=1 , MAX  00000157 

106  IF ( J. NE . JCOL (K) )  A( I , J) =A( I , J) -AI JCK*A( IROW(K) , J)  00000158 

107  CONTINUE  00000159 

C  00000160 

C  00000161 

C***********************************************************************00000162 
C*  *00000163 

C*  ORDER  SOLUTION  VALUES  (IF  ANY)  AND  CREATE  JORD  ARRAY  *00000164 

C*  *00000165 

C***********************************************************************00000166 
C  00000167 

C  00000168 

DO  108  1=1, N  00000169 

JORD ( IROW ( I ) ) =JCOL ( I )  00000170 

108  IF ( INDIC . GE . 0 )  X( JCOL(I) )=A(IROW(I) ,MAX)  00000171 

C  00000172 

C  00000173 

C***********************************************************************00000174 
C*  *00000175 

C*  ADJUST  SIGN  OF  DETERMINANT  *00000176 

C*  *00000177 

C***********************************************************************00000178 
C  00000179 

C  00000180 
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INTCH=0 

00000181 

DO  109  1=1, N-l 

00000182 

DO  109  J=I+1 , N 

00000183 

IF( JORD( J) . GE. JORD(I) )  GO  TO  109 

00000184 

JTEMP=JORD(J) 

00000185 

JORD ( J ) = JORD ( I ) 

00000186 

JORD ( I ) = JTEMP 

00000187 

INTCH=INTCH+1 

00000188 

109 

CONTINUE 

00000189 

IF ( INTCH/2*2 . NE . INTCH)  DETER=-DETER 

00000190 

C 

00000191 

C 

00000192 

C***********************************************************************00000193 

c* 

*00000194 

c* 

UNSCRAMBLE  THE  INVERSE 

*00000195 

c* 

*00000196 

C***********************************************************************00000197 

c 

00000198 

c 

00000199 

IF ( INDIC . GT . 0 )  RETURN 

00000200 

DO  111  J=1,N 

00000201 

DO  110  1=1, N 

00000202 

110 

Y ( JCOL ( I ) )=A( IROW ( I ) , J) 

00000203 

DO  111  1=1, N 

00000204 

111 

A(  I , J) =Y ( I ) 

00000205 

DO  113  1=1, N 

00000206 

DO  112  J=1 , N 

00000207 

112 

Y ( IROW( J) )=A(I, JCOL(J) ) 

00000208 

DO  113  J=1 , N 

00000209 

113 

A( I , J) =Y ( J) 

00000210 

RETURN  00000211 

C  00000212 

C  00000213 

C***********************************************************************00000214 
C*  *00000215 

C*  ERROR  MESSAGES  *00000216 

C*  *00000217 

C***********************************************************************00000218 
C  00000219 

C  00000220 

114  WRITE (UNIT=6 , FMT=201 ) N, MXNC  00000221 

RETURN  1  00000222 

C  00000223 

C  00000224 

C***********************************************************************00000225 
C*  *00000226 

C*  FORMAT  STATEMENTS  *00000227 

C*  *00000228 

C***********************************************************************00000229 
C  00000230 

C  00000231 

201  FORMAT (/,' ERROR  IN  SUBROUTINE  GJEMPS :  THE  REQUESTED  NUMBER  OF  ROWS00000232 

>,N' ,/, 'EXCEEDS  THE  DIMENSIONED  MAXIMUM  NUMBER  OF  ROWS,  MXNR  FOR' , /00000233 

>,'N  = ' , 15 , '  MXNR  =' ,15,//, 'REDIMENSION  PARAMETERS :',//, 'MXNR  =  M00000234 

>AXIMUM  NUMBER  OF  ROWS ',/,' MXNC  =  MAXIMUM  NUMBER  OF  COLUMNS ',/, 'M00000235 

>XNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  l',/,'MXNCMl  =  MAXIMUM  NUMBER00000236 

>  OF  COLUMNS  MINUS  l',//,'AND  RECOMPILE.')  00000237 

END  00000238 
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SUBROUTINE  LSPRWC(*)  00000001 
C  00000002 
C  00000003 

C***********************************************************************00000004 

C*  *00000005 
C*  LEAST  SQUARES  POLYNOMIAL  REGRESSION  WITH  CONSTRAINTS  (LSPRWC)  *00000006 
C*  *00000007 
C*  REVISION  DATE:  6  JULY  2011  *00000008 
C*  *00000009 
C***********************************************************************00000010 
C  00000011 


c* 

C*  SUBROUTINE  LSPRWC  PERFORMS  A  LEAST  SQUARES  POLYNOMIAL  REGRESSION 
C*  WITH  CONSTRAINTS;  THAT  IS,  A  SET  OF  X-Y  DATA  POINTS  IS  CURVE  FIT 
C*  WITH  AN  NP  ORDER  POLYNOMIAL  OF  THE  FORM 
C* 

C*  P(X)=B0+B1*X+B2*X**2+B3*X**3+. . . .+BNP*X**NP 

C* 

C*  WITH  ANY  POLYNOMIAL  DERIVATIVES,  ZERO  THROUGH  NP,  SPECIFIED  AT 
C*  GIVEN  X  LOCATIONS.  THE  PROCEDURE  FOLLOWED  IS  THE  METHOD  OF  LEAST 
C*  SQUARES  USING  UNDETERMINED  LAGRANGE  MULTIPLIERS. 

C* 

C*  INPUT  VARIABLES: 

C* 

C*  CV 
C*  NC 
C*  NCONST 
C*  NDATA 
C*  NP 
C*  X 
C*  XC 
C* 

C*  Y 
C* 

C*  OUTPUT  VARIABLES: 

C* 

C*  B  =  ARRAY  OF  REGRESSION  COEFFICIENTS 

C*  SDEV  =  STANDARD  DEVIATION  OF  THE  X-Y  DATA  POINTS  ABOUT  THE 
C*  REGRESSION  LINE 

C* 

C*  PARAMETERS: 

C* 

C*  MXANC  =  MAXIMUM  NUMBER  OF  ALLOWABLE  CONSTRAINTS 

C*  MXNC  =  MAXIMUM  NUMBER  OF  COLUMNS  (MXNC  =  MXNR  +  1) 

C*  MXNCM1  =  MAXIMUM  NUMBER  OF  COLUMNS  MINUS  1  (MXNCM1  =  MXNC  -  1) 

C*  MXNR  =  MAXIMUM  NUMBER  OF  ROWS 

C*  MXNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  1  (MXNRM1  =  MXNR  -  1) 

C* 

C*  REF:  CARNAHAN, B.,  LUTHER,  H.A.,  AND  WILKES,  J.O.:  APPLIED  NUMERICAL 
C*  METHODS,  NEW  YORK,  JOHN  WILEY  &  SONS,  1969,  PP.  571-584. 

C* 

C*  WE INSTOCK,  R. :  CALCULUS  OF  VARIATIONS  WITH  APPLICATIONS  TO 

C*  PHYSICS  AND  ENGINEERING,  NEW  YORK,  DOVER  PUBLICATIONS,  INC., 

C*  1974,  P.  6. 

C* 


=  ARRAY  OF  POLYNOMIAL  DERIVATIVE  VALUES 

=  ARRAY  OF  POLYNOMIAL  DERIVATIVE  ORDERS  (CONSTRAINT  ORDER) 
=  NUMBER  OF  CONSTRAINTS  OR  (XC,NC,CV)  DATA  TRIPLES 
=  NUMBER  OF  X-Y  DATA  PAIRS 
=  POLYNOMIAL  ORDER 
=  DATA  POINT  X  VALUE 

=  ARRAY  OF  POLYNOMIAL  DERIVATIVE  LOCATIONS  (CONSTRAINT  X 
VALUES) 

=  DATA  POINT  Y  VALUE 


*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 

*00000025 

*00000026 

*00000027 

*00000028 

*00000029 

*00000030 

*00000031 

*00000032 

*00000033 

*00000034 

*00000035 

*00000036 

*00000037 

*00000038 

*00000039 

*00000040 

*00000041 

*00000042 

*00000043 

*00000044 

*00000045 

*00000046 

*00000047 

*00000048 

*00000049 

*00000050 

*00000051 

*00000052 

*00000053 

*00000054 

*00000055 

*00000056 


C  00000058 

C  00000059 

IMPLICIT  REAL*8 (A-H,0-Z)  00000060 
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c 

INTEGER  COL , COLMAX , FACTRL , ROW, ROWMAX 
C 

PARAMETER  (MXANC=20 , MXNCM1=30 , MXNRM1=29) 

C 

COMMON/B/B ( 0 : MXNRM1 ) 

COMMON/C/C ( 0 : MXNRM1 , 0 : MXNCM1 ) 

COMMON /CV/CV ( MXANC ) 

COMMON / NC / NC ( MXANC ) 

COMMON/NCONST/NCONST 
COMMON / ND AT A/ ND AT A 
COMMON /NP/NP 
COMMON/ SDEV/ SDEV 
COMMON / XC / XC ( MXANC ) 

C 

c 

c* *********************************************************************  * 
C*  * 

C*  ZERO  THE  COEFFICIENT  MATRIX  AND  SOLUTION  VECTOR  * 

C*  * 

c* *********************************************************************  * 

c 

c 

ROWMAX=NP+NCONST 
COLMAX=ROWMAX+ 1 
DO  101  ROW=0, ROWMAX 
DO  101  COL=0, COLMAX 

101  C ( ROW , COL ) =0 . 0D+00 
C 

C 

c* *********************************************************************  * 
C*  * 

C*  CALCULATE  THE  COEFFICIENT  MATRIX  * 

C*  * 

c* *********************************************************************  * 

c 

c 

REWIND  1 
C 

c 

c* *********************************************************************  * 
C*  * 

C*  MINIMIZATION  EQUATIONS  -  COEFFICIENTS  OF  THE  B(J)  * 

C*  * 

c* *********************************************************************  * 

c 

c 

102  READ (UNIT=1 , END=105 , ERR=113 ) X , Y 
DO  104  JJ=0 , NP 

ROW=JJ 

DO  103  J=0 , NP 
COL=J 

103  C ( ROW , COL ) =C ( ROW , COL ) +POWER ( X , J+ J J ) 

C 

C 

c* *********************************************************************  * 
C*  * 

C*  VECTOR  OF  CONSTANTS  * 

C*  * 

c* *********************************************************************  * 


00000061 

00000062 

00000063 

00000064 

00000065 

00000066 

00000067 

00000068 

00000069 

00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 

00000088 

00000089 

00000090 

00000091 

00000092 

00000093 

00000094 

00000095 

00000096 

00000097 

00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

00000119 

00000120 
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c 

00000121 

c 

00000122 

104 

C ( ROW , COLMAX ) =C ( ROW , COLMAX ) +Y*  POWER ( X , J J ) 

00000123 

GO  TO  102 

00000124 

c 

00000125 

c 

00000126 

C***********************************************************************00000127 

c* 

*00000128 

c* 

COEFFICIENTS  OF  THE  LAMBDA ( L ) 

*00000129 

c* 

*00000130 

C***********************************************************************00000131 

c 

00000132 

c 

00000133 

105 

IF(NCONST.EQ.O)  GO  TO  108 

00000134 

DO  106  JJ=0 , NP 

00000135 

ROW=JJ 

00000136 

F1=FACTRL( JJ) 

00000137 

DO  106  L=1 , NCONST 

00000138 

IF( JJ.LT.NC(L) )  GO  TO  106 

00000139 

COL=NP+L 

00000140 

F2=F1/FACTRL ( JJ-NC (L) ) /2 . 0D+00 

00000141 

C (ROW, COL) =F2*POWER(XC (L) , JJ-NC (L) ) 

00000142 

106 

CONTINUE 

00000143 

C 

00000144 

C 

00000145 

C***********************************************************************00000146 

c* 

*00000147 

c* 

CONSTRAINT  EQUATIONS  -  COEFFICIENTS  OF  THE  B(J) 

*00000148 

c* 

*00000149 

C***********************************************************************00000150 

c 

00000151 

c 

00000152 

DO  107  L=l, NCONST 

00000153 

ROW=NP+L 

00000154 

C ( ROW , COLMAX ) =CV ( L ) 

00000155 

DO  107  J=NC (L) , NP 

00000156 

COL=J 

00000157 

F1=FACTRL( J) /FACTRL( J-NC(L) ) 

00000158 

107 

C (ROW, COL) =F1*P0WER(XC (L) , J-NC (L) ) 

00000159 

C 

00000160 

C 

00000161 

C***********************************************************************00000162 

c* 

*00000163 

c* 

COMPUTE  THE  SOLUTION  VECTOR 

*00000164 

c* 

*00000165 

C***********************************************************************00000166 

c 

00000167 

c 

00000168 

108 

CALL  GJEMPS ( 1 . 0D-20 , +1 , ROWMAX+1 , DETER, * 112 ) 

00000169 

IF ( DETER . EQ . 0 . 0D+00 )  GO  TO  114 

00000170 

C 

00000171 

C 

00000172 

C***********************************************************************00000173 

c* 

*00000174 

c* 

COMPUTE  THE  STANDARD  DEVIATION 

*00000175 

c* 

*00000176 

C***********************************************************************00000177 

c 

00000178 

c 

00000179 

SDEV=0 . 0D+00 

00000180 
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IDENOM=NDATA- (NP+1 ) -NCONST-1 

00000181 

IF ( IDENOM. LE . 0 )  RETURN 

00000182 

DENOM=IDENOM 

00000183 

REWIND  1 

00000184 

109 

READ (UNIT=1 , END=111 , ERR=113 ) X,  Y 

00000185 

PX=0 . 0D+00 

00000186 

DO  110  J=NP ,1,-1 

00000187 

110 

PX=(PX+B(J) )*X 

00000188 

PX=PX+B(0) 

00000189 

SDEV=SDEV+ ( Y-PX ) *  *  2 

00000190 

GO  TO  109 

00000191 

111 

SDEV=DSQRT ( SDEV/DENOM) 

00000192 

RETURN 

00000193 

C 

00000194 

C 

00000195 

C***********************************************************************00000196 

c* 

*00000197 

c* 

ERROR  MESSAGES 

*00000198 

c* 

*00000199 

C***********************************************************************00000200 

c 

00000201 

c 

00000202 

112 

WRITE (UNIT=6 , FMT=201 ) 

00000203 

RETURN  1 

00000204 

113 

WRITE (UNIT=6 , FMT=202 ) 

00000205 

RETURN  1 

00000206 

114 

WRITE (UNIT=6 , FMT=203 ) 

00000207 

RETURN  1 

00000208 

C 

00000209 

c 

00000210 

C***********************************************************************00000211 

c* 

*00000212 

c* 

FORMAT  STATEMENTS 

*00000213 

c* 

*00000214 

c 

c 

201 

202 

203 


FORMAT ( ' SUBROUTINE 
FORMAT ( / ,  ' ERROR  IN 
>IT  =  1' ) 

FORMAT ( ' SUBROUTINE 
END 


GJEMPS  WAS  CALLED 
SUBROUTINE  LSPRWC 

GJEMPS  WAS  CALLED 


FROM  SUBROUTINE  LSPRWC') 
UNFORMATTED  READ  ERROR 

FROM  SUBROUTINE  LSPRWC') 


00000216 
00000217 
00000218 
ON  UN00000219 
00000220 
00000221 
00000222 
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non  h  o  o 


SUBROUTINE  PES(*) 
C 
C 


C*  * 

C*  POLYNOMIAL  EVALUATION  SUBROUTINE  (PES)  * 

C*  * 

C*  REVISION  DATE:  3  MAY  2002  * 

C*  * 


c 


c* *********************************************************************  * 

c* 

C*  GIVEN  X  AND  THE  COEFFICIENTS  B(0) ,B(1) , . . . . ,B(NP) ,  SUBROUTINE  PES 
C*  EVALUATES  THE  K-TH  DERIVATIVE  OF  POLYNOMIALS  OF  THE  FORM: 

C* 

C*  P (X) =B (0) +B ( 1 ) *X+B (2 ) *X**2+ . . . .+B(NP) *X**NP 

C* 

C*  INPUT  VARIABLES: 

C* 

C*  B  =  ARRAY  OF  POLYNOMIAL  COEFFICIENTS 

C*  NP  =  POLYNOMIAL  ORDER 

C*  X  =  X  VALUE  FOR  EVALUATION  OF  THE  POLYNOMIAL 

C* 

C*  OUTPUT  VARIABLES: 

C* 


c* 

R 

=  VALUE 

AT 

X 

OF 

THE 

c* 

S 

=  VALUE 

AT 

X 

OF 

THE 

c* 

POLYNOMIAL 

c* 

YK 

=  VALUE 

AT 

X 

OF 

THE 

C* 

C*  PARAMETERS: 

C* 

C*  MXNR  =  MAXIMUM  NUMBER  OF  ROWS 

C*  MXNRM1  =  MAXIMUM  NUMBER  OF  ROWS  MINUS  1  (MXNRM1  =  MXNR 
C* 

C*  *********************************************************************  * 


1) 


c 

c 

c 

c 

c 

c 


IMPLICIT  REAL*8 (A-H,0-Z) 

INTEGER  FACTRL 

PARAMETER  (MXNRM1=29) 

COMMON/B/B ( 0 : MXNRM1 ) 
COMMON /NP/NP 

DIMENSION  YK ( 0 : 9 ) 


REWIND  2 
REWIND  3 
S=0 . 0E+00 
N=0 


READ (UNIT=2 , END=106 , ERR=107 ) X 
N=N+1 


*********************************************************************** 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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C*  * 

C*  EVALUATE  THE  NK-TH  DERIVATIVE  OF  THE  POLYNOMIAL  (0  TO  NP)  * 

C*  * 

c* *********************************************************************  * 

c 

c 

DO  103  NK=0,NP 
YK(NK)=0.0D+00 
DO  102  J=NP, NK+1 , -1 

102  YK(NK) = ( YK(NK) +FACTRL ( J) *B ( J) /FACTRL ( J-NK) ) *X 

103  YK(NK)=YK(NK)+FACTRL(NK) *B(NK) 

C 

C 

C*  *********************************************************************  * 

C*  * 

C*  EVALUATE  THE  RADIUS  OF  CURVATURE  OF  THE  POLYNOMIAL  * 

C*  * 

c* *********************************************************************  * 

c 

c 

R=1 . OE+32 

IF (NP . LT . 2 )  GO  TO  104 

IF ( YK ( 2 ) . EQ . 0 . 0E+00 )  GO  TO  104 

R=( (1.0E+00+YK(l)*YK(l) ) **1.5) /DABS (YK(2) ) 

104  CONTINUE 
C 

C 

C*  *********************************************************************  * 

C*  * 

C*  CALCULATE  THE  ARC  LENGTH  * 

C*  * 

c* *********************************************************************  * 

c 

c 

c 

IF(N.EQ.l)  GO  TO  105 
DeltaX=X-Xsave 
DeltaY=YK ( 0 ) -Ysave 

DeltaS=DSQRT(DeltaX*DeltaX+DeltaY*DeltaY) 

S=S+DeltaS 

105  Xsave=X 
Ysave=YK(0) 

WRITE (UNIT=3 , ERR=108 ) X, (YK(K) ,K=0,NP) ,R,S 
GO  TO  101 
C 
C 

106  END  FILE  3 
RETURN 

C 

C 

C*  *********************************************************************  * 

C*  * 

C*  ERROR  MESSAGES  * 

C*  * 


c 

c 

107  WRITE (UNIT=6,FMT=201) 
RETURN  1 

108  WRITE (UNIT=6 , FMT=202 ) 


00000061 

00000062 

00000063 

00000064 

00000065 

00000066 

00000067 

00000068 

00000069 

00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 

00000088 

00000089 

00000090 

00000091 

00000092 

00000093 

00000094 

00000095 

00000096 

00000097 

00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

00000119 

00000120 
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RETURN  1  00000121 
C  00000122 
C  00000123 

C***********************************************************************00000124 

C*  *00000125 
C*  FORMAT  STATEMENTS  *00000126 
C*  *00000127 

C***********************************************************************00000128 


c 

c 


201 

FORMAT ( / 

>=  2') 

202 

FORMAT ( / 

>  =  3') 

END 

'ERROR  IN  SUBROUTINE  PES: 
'ERROR  IN  SUBROUTINE  PES: 


00000129 

00000130 

UNFORMATTED  READ  ERROR  ON  UNIT  00000131 

00000132 

UNFORMATTED  WRITE  ERROR  ON  UNIT00000133 

00000134 

00000135 
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SUBROUTINE  PODS ( * ) 

C 

C 

C*  ********************************************************************  * 

c* 

C*  POLYNOMIAL  ORDER  DEFINITION  SUBROUTINE  (PODS) 

C* 

C*  REVISION  DATE:  11  SEPTEMBER  2007 

C* 

C*  ********************************************************************  * 

c 

c* ********************************************************************  * 

c* 

C*  SUBROUTINE  PODS  DEFINES  THE  POLYNOMIAL  ORDER  FOR  PROGRAM  LSPRWC. 

C* 

C*  INPUT  VARIABLES: 

C* 

C*  NDATA  =  NUMBER  OF  X,Y  DATA  POINTS  TO  BE  FIT 
C*  NCONST  =  NUMBER  OF  CONSTRAINTS 
C* 

C*  INPUT/OUTPUT  VARIABLES: 

C* 

C*  NP  =  POLYNOMIAL  ORDER 

C* 

C*  ********************************************************************  * 

c 

c* ********************************************************************  * 

c* 

C*  PARAMETERS: 

C* 

C*  MXNAC  =  MAXIMUM  NUMBER  OF  ALLOWABLE  CONSTRAINTS 
C* 

C*  ********************************************************************  * 

c 

c 

IMPLICIT  REAL*8 (A-H,0-Z) 

C 

PARAMETER  (MXNAC=20 , MXNR=30 , MXNRM1=29 ) 

C 

COMMON / NC / NC ( MXNAC ) 

COMMON/NCONST/NCONST 
COMMON / NDATA/ NDATA 
COMMON /NP/NP 
C 

DATA  NP/3/ 

C 

NAMEL I ST / PARM/NP 
C 
C 

c* ********************************************************************  * 

c* 

C*  INPUT  THE  POLYNOMIAL  ORDER 

C* 

C*  ********************************************************************  * 

c 

c 

NPMAX=MINO ( 9 , NDATA-NCONST-1 ) 

101  WRITE (UNIT=6,FMT=201)NPMAX 
CALL  NDWUS ( ' PARM ' , 6 , * 103 ) 

CALL  NDWUSI( 'NP' ,1,1,1, NP, *103) 


00000001 

00000002 

00000003 

*00000004 

*00000005 

*00000006 

*00000007 

*00000008 

*00000009 

*00000010 

00000011 

*00000012 

*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 

*00000025 

00000026 

*00000027 

*00000028 

*00000029 

*00000030 

*00000031 

*00000032 

*00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

*00000050 

*00000051 

*00000052 

*00000053 

*00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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102 

C 

C 


CALL  NDWUSE ( *103 ) 

CALL  NDRUS ( ' PARM' ,4, *101, *104) 

READ (UNIT=4 , NML=PARM, ERR=101 ) 

CALL  NDRUSE(*101, *104) 

IF (NP . LT . 0 . OR. NP . GT . 9 )  GO  TO  105 
IF (NP . GT . NPMAX)  GO  TO  106 
DO  102  1=1 , NCONST 

IF (NC ( I ) . LT . 0 . OR. NC ( I ) . GT . NP)  GO  TO  107 

CONTINUE 

RETURN 


00000061 

00000062 

00000063 

00000064 

00000065 

00000066 

00000067 

00000068 

00000069 

00000070 

00000071 


00000072 

C***********************************************************************00000073 

C*  *00000074 

C*  ERROR  MESSAGES  *00000075 

C*  *00000076 

C***********************************************************************00000077 


c 

c 

103 

104 

105 

106 

107 

108 


WRITE (UNIT=6 , FMT=204 ) 

GO  TO  108 

WRITE (UNIT=6 , FMT=205 ) 

GO  TO  108 

WRITE (UNIT=6 , FMT=206 ) NP 
GO  TO  108 

WRITE (UNIT=6 , FMT=207 ) NP , NCONST , NDATA, NPMAX 
GO  TO  108 

WRITE (UNIT=6 , FMT=208 ) NC ( I ) , I , NP 
WRITE (UNIT=6 , FMT=202 ) 

READ (UNIT=5 , FMT=203 ) CCV 
RETURN  1 


00000078 
00000079 
00000080 
00000081 
00000082 
00000083 
00000084 
00000085 
00000086 
00000087 
00000088 
00000089 
00000090 
00000091 

C  00000092 

C  00000093 

C***********************************************************************00000094 

C*  *00000095 

C*  FORMAT  STATEMENTS  *00000096 

C*  *00000097 

C***********************************************************************00000098 

C  00000099 

C  00000100 

201  FORMAT (/,' INPUT  THE  POLYNOMIAL  ORDER  BY  NAMELIST  WHERE NP  00000101 

>  =  THE  POLYNOMIAL  ORDER' ,11,' NOTE :  NP  MUST  BE  IN  THE  RANGE  OF  0  T000000102 

>  9. ' ,/,T7, 'THERE  MUST  BE  AT  LEAST  NP+NCONST+1  X-Y  DATA  POINTS 00000103 
>T7 , ' THE  MAXIMUM  VALUE  FOR  NP  IS  ', 15 CURRENT  VALUES  ARE :') 00000104 

FORMAT (/, T19, '-  ENTER/RETURN  TO  CONTINUE  -')  00000105 

FORMAT (Al)  00000106 

FORMAT ( 'SUBROUTINE  NDWUS  WAS  CALLED  FROM  SUBROUTINE  PODS')  00000107 

FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  SUBROUTINE  PODS')  00000108 

FORMAT (/,' ERROR  IN  SUBROUTINE  PODS:  NP  =  ',15,'  IS  NOT  ALLOWED .', /00000109 
>, ' THE  POLYNOMIAL  ORDER  MUST  BE  GREATER  THAN  OR  EQUAL  TO  0  AND  LESS00000110 

>  THAN  OR  EQUAL',/, 'TO  9.')  00000111 

FORMAT (/,' ERROR  IN  SUBROUTINE  PODS:  NP  =  ',15,'  IS  NOT  ALLOWED  FOROOOOOH2 

>  NCONST  =  ',15,'  AND ',/,' NDATA  =  ',15,'.  THERE  MUST  BE  AT  LEAST  N00000113 

>P+2*NCONST  X-Y  DATA  POINTS.  THE ',/, 'MAXIMUM  VALUE  FOR  NP  IS  ',15,00000114 

>'.')  00000115 

FORMAT (/,' ERROR  IN  SUBROUTINE  PODS:  NC(',I3,')  =  ',15,'  IS  NOT  ALL00000116 
>OWED  FOR  NP  =  ' ,15, '.',/, 'THE  CONSTRAINT  ORDER  MUST  BE  GREATER  THA00000117 
>N  OR  EQUAL  TO  0  AND  LESS  THAN  OR  EQUAL  TO  NP.')  00000118 

END  00000119 


202 

203 

204 

205 

206 


207 


208 
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FUNCTION  POWER (X,N) 


00000001 

C  00000002 

C  00000003 

C***********************************************************************00000004 

C*  *00000005 

C*  INTEGER  POWERS  OF  REAL  NUMBERS  FUNCTION  SUBPROGRAM  (POWER)  *00000006 

C*  *00000007 

C*  REVISION  DATE:  15  JUNE  1989  *00000008 

C*  *00000009 

C***********************************************************************00000010 


C  00000011 

C  00000012 

REAL *8  POWER, X  00000013 

POWER=l . 0D+00  00000014 

IF(N.EQ.O)  RETURN  00000015 

DO  101  1=1, N  00000016 

101  POWER=POWER*X  00000017 

END  00000018 
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SUBROUTINE  QPS ( * ) 

C 

C 

C*  ********************************************************************  * 

c* 

C*  QUICK  PLOT  SUBROUTINE  (QPS) 

C* 

C*  REVISION  DATE:  1  JULY  2011 

C* 

C*  ********************************************************************  * 

c 

c* ********************************************************************  * 

c* 

C*  SUBROUTINE  QPS  WRITES  DATA  FILES  TO  CONSTRUCT  A  QUICK  PLOT  OF 
C*  RESULTS  FOR  PROGRAM  LSPRWC  USING  THE  PostScript  CARTESIAN 
C*  COORDINATE  PLOT  PROGRAM  PSCCPP. 

C* 

C*  INPUT  VARIABLES: 

C* 

C*  NCONST  =  NUMBER  OF  CONSTRAINTS 

C*  X  =  X  VALUE  FOR  EVALUATION  OF  THE  POLYNOMIAL 

C*  XC  =  ARRAY  OF  CONSTRAINT  X  VALUES 

C*  Y  =  Y  VALUE  FROM  EVALUATION  OF  THE  POLYNOMIAL  AT  X 

C* 

C*  PARAMETERS: 

C* 

C*  MXNAC  =  MAXIMUM  NUMBER  OF  ALLOWABLE  CONSTRAINTS 
C* 

C*  ********************************************************************  * 

c 

c 

IMPLICIT  REAL*8 (A-H,0-Z) 

C 

CHARACTER* 10  File 
C 

PARAMETER  (MXNAC=20 , One=l . 0E+00 , Zero=0 . 0E+00 ) 

C 

COMMON/NCONST/NCONST 
COMMON / XC / XC ( MXNAC ) 

C 

DIMENSION  File(3) 

C 

DATA  File/'@lsprwc.cl' , '@lsprwc.ed' , ' @lsprwc .dp ' / , Norm/0/ 

C 

NAMELIST/PARM/Norm, Xmax , Xmin , Ymax , Ymin 
C 

FNX (Q)=( Q-Xmin ) / ( Xmax-Xmin ) 

FNY (Q)=( Q- Ymin ) / ( Ymax- Ymin ) 

C 

C 


00000001 

00000002 

00000003 

*00000004 

*00000005 

*00000006 

*00000007 

*00000008 

*00000009 

*00000010 

00000011 

*00000012 

*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 

*00000025 

*00000026 

*00000027 

*00000028 

*00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 


C* 

C* 

C* 


FIND  THE  RANGE  OF  THE  DATA 


*00000052 

*00000053 

*00000054 


C  00000056 

C  00000057 

Xmin=+1.0D+75  00000058 

Xmax=-1.0D+75  00000059 

Ymin=+1.0D+75  00000060 
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Ymax=-1.0D+75 
REWIND  1 

101  READ (UNIT=1 , END=102 , ERR=116 ) X , Y 
Xmin=DMINl ( X , Xmin ) 

Xmax=DMAX 1  ( X , Xmax ) 

Ymin=DMINl ( Y , Ymin ) 

Ymax=DMAXl ( Y , Ymax ) 

GO  TO  101 

102  REWIND  3 

103  READ (UNIT=3 , END=104 , ERR=117 ) X , Y 
Xmin=DMINl ( X , Xmin ) 

Xmax=DMAX 1  ( X , Xmax ) 

Ymin=DMINl ( Y , Ymin ) 

Ymax=DMAX 1  ( Y , Ymax ) 

GO  TO  103 

104  IF ( NCONST . EQ . 0 )  GO  TO  106 
DO  105  1=1, NCONST 
Xmin=DMINl (XC ( I ) , Xmin ) 

105  Xmax=DMAX 1 ( XC ( I ) , Xmax ) 

106  IF ( Xmin. GE. Xmax)  RETURN 
IF ( Ymin. GE. Ymax)  RETURN 

107  WRITE (UNIT=6 , FMT=201 ) 

CALL  NDWUS ( ' PARM ' , 6 , * 1 1 8 ) 

CALL  NDWUSI( 'Norm' , 1 , 1 , 1 , Norm, *118 ) 
CALL  NDWUSD (' Xmax ' , 1, 1, 1, Xmax, *118) 
CALL  NDWUSD (' Xmin ' , 1 , 1 , 1 , Xmin, *118 ) 
CALL  NDWUSD (' Ymax ' , 1 , 1 , 1 , Ymax, *118 ) 
CALL  NDWUSD (' Ymin ' , 1 , 1 , 1 , Ymin, *118 ) 
CALL  NDWUSE(*118) 

CALL  NDRUS ( ' PARM' ,4, *107, *119) 

READ (UNIT=4 , NML=PARM, ERR=107 ) 

CALL  NDRUSE(*107 , *119) 

C 

C 

c** 

c* 
c* 
c* 

c** 

c 
c 


****************************************************************** 


MARK  THE  CONSTRAINT  LOCATIONS 


108 

109 

110 
C 

C 

C** 

c* 


****************************************************************** 


N=1 

OPEN (UNIT=4 , FILE=File (N) , STATUS= ' UNKNOWN ' , FORM= ' FORMATTED ' , 
>ERR=120) 

IF ( NCONST. EQ.0)  GO  TO  109 
DO  108  1=1, NCONST 

IF ( Norm . EQ . 0 )  WRITE (UNIT=4 , FMT=202 , ERR=12 1 ) XC ( I ) , Ymin , INT (Zero ) 
IF (Norm. EQ . 1 )  WRITE (UNIT=4 , FMT=202 , ERR=121 ) FNX (XC ( I ) ) , Zero , 
>INT(Zero) 

IF ( Norm . EQ . 0 )  WRITE (UNIT=4 , FMT=202 , ERR=12 1 ) XC ( I ) , Ymax , INT ( One ) 
IF ( Norm . EQ . 1 )  WRITE (UNIT=4 , FMT=202 , ERR=12 1 ) FNX ( XC ( I ) ) , One , 
>INT(One) 

CONTINUE 
GO  TO  110 

WRITE (UNIT=4 , FMT=202 , ERR=121 ) Zero , Zero , INT ( Zero) 

CLOSE (UNIT=4 , STATUS= ' KEEP ' , ERR=122 ) 


****************************************************************** 


00000061 

00000062 

00000063 

00000064 

00000065 

00000066 

00000067 

00000068 

00000069 

00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 

00000088 

00000089 

00000090 

00000091 

00000092 

00000093 

00000094 

***00000095 

*00000096 

*00000097 

*00000098 

***00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

***00000119 

*00000120 
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C*  WRITE  THE  EVALUATED  DATA  * 

C*  * 

c* *********************************************************************  * 

c 

c 

N=2 

OPEN (UNIT=4 , FILE=File (N) , STATUS= ' UNKNOWN ' , FORM= ' FORMATTED ' , 
>ERR=120) 

111  REWIND  3 

112  READ (UNIT=3 , END=113 , ERR=117 ) X, Y 

IF (Norm. EQ . 0 )  WRITE (UNIT=4 , FMT=202 , ERR=121 ) X , Y 
IF(Norm.EQ.l)  WRITE (UNIT=4 , FMT=202 , ERR=121 ) FNX (X) , FNY ( Y) 

GO  TO  112 

113  CLOSE (UNIT=4,STATUS=' KEEP' ,ERR=122) 

C 

C 

c* *********************************************************************  * 
C*  * 

C*  MARK  THE  X-Y  DATA  POINTS  TO  BE  FITTED  * 

C*  * 

c* *********************************************************************  * 

c 

c 

N=3 

OPEN (UNIT=4 , FILE=File (N) , STATUS= ' UNKNOWN ' , FORM= ' FORMATTED ' , 
>ERR=120) 

REWIND  1 

114  READ (UNIT=1 , END=115 , ERR=116 )X, Y 

IF ( Norm . EQ . 0 )  WRITE (UNIT=4 , FMT=202 , ERR=12 1 ) X , Y 
IF(Norm.EQ.l)  WRITE (UNIT=4 , FMT=202 , ERR=121 ) FNX (X) , FNY ( Y) 

GO  TO  114 

115  CLOSE (UNIT=4,STATUS=' KEEP' ,ERR=122) 

C 

C 

c* *********************************************************************  * 
C*  * 

C*  CONTSTRUCT  THE  PLOT  * 

C*  * 

c* *********************************************************************  * 

c 

c 

RETURN 

C 

C 

c* *********************************************************************  * 
C*  * 

C*  ERROR  MESSAGES  * 

C*  * 

c* *********************************************************************  * 

c 

c 

116  WRITE (UNIT=6 , FMT=205 ) 

RETURN  1 

117  WRITE (UNIT=6 , FMT=206 ) 

RETURN  1 

118  WRITE (UNIT=6 , FMT=207 ) 

RETURN  1 

119  WRITE (UNIT=6 , FMT=208 ) 

RETURN  1 

120  WRITE (UNIT=6,FMT=209) File (N) 


00000121 

00000122 

00000123 

00000124 

00000125 

00000126 

00000127 

00000128 

00000129 

00000130 

00000131 

00000132 

00000133 

00000134 

00000135 

00000136 

00000137 

00000138 

00000139 

00000140 

00000141 

00000142 

00000143 

00000144 

00000145 

00000146 

00000147 

00000148 

00000149 

00000150 

00000151 

00000152 

00000153 

00000154 

00000155 

00000156 

00000157 

00000158 

00000159 

00000160 

00000161 

00000162 

00000163 

00000164 

00000165 

00000166 

00000167 

00000168 

00000169 

00000170 

00000171 

00000172 

00000173 

00000174 

00000175 

00000176 

00000177 

00000178 

00000179 

00000180 


B-32 


RETURN  1 

121  WRITE (UNIT=6 , FMT=210) File (N) 
RETURN  1 

122  WRITE (UNIT=6,FMT=211) File (N) 
RETURN  1 

C 

C 


00000181 

00000182 

00000183 

00000184 

00000185 

00000186 

00000187 


C***********************************************************************00000188 

C*  *00000189 

C*  FORMAT  STATEMENTS  *00000190 

C*  *00000191 

C***********************************************************************00000192 

C  00000193 

C  00000194 

201  FORMAT (/,' INPUT  THE  PLOT  VARIABLES  BY  NAMELIST  WHERE Norm  00000195 

>=  RENormILIZATION  CONTROL T8 ,' =  0,  TO  NOT  NormALIZE  THE  PLOT  DA00000196 
>TA',/,T8,'=  1,  TO  NormALIZE  THE  PLOT  DATA ' , / , ' Xmax  =  MAXIMUM  X-V00000197 
>ALUE ' , / / ' Xmin  =  MINIMUM  X-VALUE '  ,  / ,  ' Ymax  =  MAXIMUM  Y-VALUE '  ,  / , 00000198 
>' Ymin  =  MINIMUM  Y-VALUE ' ,11,' CURRENT  VALUES  ARE:')  00000199 

202  FORMAT (2E 15 .6,115)  00000200 

203  FORMAT ( 'lsprwc.l' )  00000201 

204  FORMAT ( 'lsprwc.ps' )  00000202 

205  FORMAT (/,' ERROR  IN  SUBROUTINE  QPS :  UNFORMATTED  READ  ERROR  ON  UNIT100000203 

>')  00000204 

206  FORMAT (/,' ERROR  IN  SUBROUTINE  QPS:  UNFORMATTED  READ  ERROR  ON  UNIT300000205 

>')  00000206 

207  FORMAT ( 'SUBROUTINE  NDWUS  WAS  CALLED  FROM  SUBROUTINE  QPS')  00000207 

208  FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  SUBROUTINE  QPS')  00000208 

209  FORMAT (/,' ERROR  IN  SUBROUTINE  QPS:  OPEN  ERROR  ON  UNIT  4,  STATUS  =  00000209 

>" UNKNOWN" ,  FILE  =',A10)  00000210 

210  FORMAT (/,' ERROR  IN  SUBROUTINE  QPS:  WRITE  ERROR  ON  UNIT  4,  FILE  =  '00000211 

>, A10)  00000212 

211  FORMAT (/,' ERROR  IN  SUBROUTINE  QPS:  CLOSE  ERROR  ON  UNIT  4,  FILE  =',00000213 

>A10)  00000214 

END  00000215 
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c 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c 


SUBROUTINE  XDPEDS (LCV2 , * , * ) 


********************************************************************** 

* 

X  DATA  POINTS  FOR  EVALUATION  DEFINITION  SUBROUTINE  (XDPEDS)  * 

* 

REVISION  DATE:  1  JULY  2011  * 

* 

********************************************************************** 

********************************************************************** 


SUBROUTINE  XDPEDS  DEFINES  THE  X  COORDINATE  SET  FOR  EVALUATION  OF  * 
THE  RESULTANT  LEAST  SQUARES  POLYNOMIAL  SET.  * 

* 

INPUT/OUTPUT  VARIABLES:  * 

* 

X  =  ARRAY  OF  X  COORDINATES  FOR  EVALUATION  OF  THE  LEAST  SQUARES  * 

POLYNOMIAL  SET  * 

* 

OUTPUT  VARIABLES :  * 

* 

NEVAL  =  NUMBER  OF  X  COORDINATES  FOR  EVALUATION  * 

* 

PARAMETERS :  * 

* 

MNXCE  =  MAXIMUM  NUMBER  OF  X  COORDINATES  FOR  POLYNOMIAL  EVALUATION  * 

* 

********************************************************************** 


IMPLICIT  REAL*8 (A-H,0-Z) 

CH ARACTE R  CCV*l,DFmt*3,File*80,Fmt*80 

PARAMETER  (R_NaN=-999 .E+00) 

PARAMETER  (MNXCE=101) 

COMMON / NEVAL / NEVAL 
COMMON/XINT/XINT 
COMMON / XMAX / XMAX 
COMMON/XMIN/XMIN 

DIMENSION  X (MNXCE) 

DATA  DFmt/ ' ( * ) ' / ,NEVAL/0/ 

NAMELIST/PARM/X, XINT , XMAX , XMIN 


********************************************************************** 

* 

INPUT  THE  X  DATA  POINTS  FOR  EVALUATION  * 

OF  THE  LEAST  SQUARES  POLYNOMIAL  * 

* 

********************************************************************** 


IF ( NEVAL . EQ . 0 )  WRITE (UNIT=6 , FMT=201 ) 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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IF ( NEVAL . NE . 0 )  WRITE (UNIT=6 , FMT=202 )  00000061 

READ (UNIT=5 , FMT=203 ) LCV1  00000062 

REWIND  2  00000063 

LCV2=1  00000064 

GO  TO  (101,107,112,116) ,LCV1  00000065 

C  00000066 

C  00000067 


C***********************************************************************00000068 


C*  *00000069 
C*  NAMELIST  INPUT  *00000070 
C*  *00000071 

C***********************************************************************00000072 


c 

c 

101 

102 

103 


104 

105 

106 

C 

C 


DO  102  1=1 ,MNXCE 
X ( I ) =R_NaN 

WRITE (UNIT=6 , FMT=204 )MNXCE 
WRITE (UNIT=6 , FMT=205 ) 

CALL  NDRUS ( ' PARM' , 4 , *103 , *117 ) 
READ (UNIT=4 , NML=PARM, ERR=103 ) 
CALL  NDRUSE ( *103 , *117 ) 

DO  104  NE VAL=MNXCE ,1,-1 

IF (X (NEVAL) .NE.RNaN)  GO  TO  105 

CONTINUE 

DO  106  1=1, NEVAL 

WRITE (UNIT=2 , ERR=122 ) X ( I ) 

GO  TO  115 


00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 

00000088 


00000089 

C***********************************************************************00000090 

C*  *00000091 

C*  FORMATTED  DISK  FILE  INPUT  *00000092 

C*  *00000093 

C***********************************************************************00000094 


c 

c 

107 


108 


109 


110 


WRITE (UNIT=6 , FMT=206 ) 

WRITE (UNIT=6 , FMT=207 ) 

READ (UNIT=5 , FMT=208 ) File 
WRITE (UNIT=6 , FMT=209 ) 

WRITE (UNIT=6 , FMT=210) 

CALL  DFCUS ( FMT , DFmt , * 108 ) 

OPEN (UNIT=4 , FILE=File, STATUS= ' OLD ' , ERR=118 ) 

REWIND  4 

NEVAL=0 

IF ( Fmt (1:3) . EQ . DFmt )  GO  TO  110 

READ ( UNIT=4 , FMT=Fmt , END=1 1 1 , ERR=1 1 9 ) X ( 1 ) 

WRITE (UNIT=2 , ERR=122 ) X ( 1 ) 

NE VAL=NE VAL+ 1 
GO  TO  109 

READ (UNIT=4 , FMT=* , END=1 1 1 , ERR=1 19 ) X ( 1 ) 

WRITE (UNIT=2 , ERR=122 ) X ( 1 ) 

NE VAL=NE VAL+ 1 
GO  TO  110 

CLOSE (UNIT=4 , STATUS= ' KEEP ' , ERR=120 ) 

GO  TO  115 


111 

C 
C 

C*  *00000120 


00000095 

00000096 

00000097 

00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 
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C*  MINIMUM,  MAXIMUM,  AND  INTERVAL  VALUES  *00000121 

C*  *00000122 

C***********************************************************************00000123 
C  00000124 

C  00000125 

112  WRITE (UNIT=6,FMT=211)  00000126 

CALL  NDWUS( 'PARM' ,6, *121)  00000127 

CALL  NDWUSD( 'XINT' , 1,1,1, XINT, *121)  00000128 

CALL  NDWUSD ( ' XMAX ',1,1,1, XMAX ,*121)  00000129 

CALL  NDWUSD( 'XMIN' , 1,1,1, XMIN, *121)  00000130 

CALL  NDWUSE(*121)  00000131 

CALL  NDRUS ( ' PARM' , 4 , *112 , *117 )  00000132 

READ (UNIT=4 , NML=PARM, ERR=112 )  00000133 

CALL  NDRUSE (*112, *117)  00000134 

NEVAL=0  00000135 

X(1)=XMIN  00000136 

DUMl=XMAX+XINT/2  00000137 

113  IF(X(1) . GT.DUM1)  GO  TO  114  00000138 

WRITE (UNIT=2 , ERR=122 ) X ( 1 )  00000139 

NE VAL=NE VAL+ 1  00000140 

DUM3  =X ( 1 )  00000141 

X ( 1 ) =XMIN+NEVAL*XINT  00000142 

GO  TO  113  00000143 

114  IF(DUM3 . EQ.XMAX)  GO  TO  115  00000144 

WRITE (UNIT=2 , ERR=122 ) XMAX  00000145 

NE VAL=NE VAL+ 1  00000146 

115  END  FILE  2  00000147 

116  RETURN  00000148 

C  00000149 

C  00000150 

C***********************************************************************00000151 
C*  *00000152 

C*  ERROR  MESSAGES  *00000153 

C*  *00000154 

C***********************************************************************00000155 
C  00000156 

C  00000157 

117  WRITE (UNIT=6 , FMT=212 )  00000158 

GO  TO  124  00000159 

118  WRITE (UNIT=6 , FMT=213 )  00000160 

CALL  CBUS06(6, File, *123)  00000161 

RETURN  1  00000162 

119  WRITE (UNIT=6 , FMT=214 )  00000163 

CALL  CBUS06(6, File, *123)  00000164 

RETURN  1  00000165 

120  WRITE (UNIT=6 , FMT=215 )  00000166 

CALL  CBUS06(6, File, *123)  00000167 

RETURN  1  00000168 

121  WRITE (UNIT=6 , FMT=216 )  00000169 

GO  TO  124  00000170 

122  WRITE (UNIT=6 , FMT=217 )  00000171 

RETURN  1  00000172 

123  WRITE (UNIT=6 , FMT=218 )  00000173 

RETURN  1  00000174 

124  WRITE (UNIT=6 , FMT=219 )  00000175 

READ (UNI T= 5 , FMT=220 ) CCV  00000176 

RETURN  2  00000177 

C  00000178 

C  00000179 

C***********************************************************************00000180 


117  WRITE (UNIT=6 , FMT=212 ) 

GO  TO  124 

118  WRITE (UNIT=6 , FMT=213 ) 
CALL  CBUS06 (6, File, *123) 
RETURN  1 

119  WRITE (UNIT=6 , FMT=214 ) 
CALL  CBUS06 (6 , File, *123) 
RETURN  1 

120  WRITE (UNIT=6 , FMT=215 ) 
CALL  CBUS06 (6, File, *123) 
RETURN  1 

121  WRITE (UNIT=6,FMT=2 16) 

GO  TO  124 

122  WRITE (UNIT=6 , FMT=217 ) 
RETURN  1 

123  WRITE (UNIT=6,FMT=2 18) 
RETURN  1 

124  WRITE (UNIT=6 , FMT=219 ) 
READ (UNI T= 5 , FMT=220 ) CCV 
RETURN  2 
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c* 

c* 

c* 

c*  *  *  * 

c 

c 

201 


202 


203 

204 

205 

206 

207 

208 

209 

210 

211 


212 

213 

214 

215 

216 

217 

218 

219 

220 


FORMAT  STATEMENTS 


*00000181 
*00000182 
*00000183 

*******************************************************************00000184 

00000185 

00000186 

FORMAT (/,' SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X  DATA  TO  BE  EVALUATE00000187 
>D  FROM  THE  FOLLOWING' ,/, 'LIST: 1,  FOR  KEYBOARD  INPUT  VIA  NAME00000188 
>LIST ' , / , ' 2 ,  FOR  FORMATTED  DISK  FILE  INPUT',/, '3,  FOR  INPUT  OF  MINI00000189 
>MUM,  MAXIMUM,  AND  INTERVAL  VALUES')  00000190 

FORMAT (/,' SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X  DATA  TO  BE  EVALUATE00000191 
>D  FROM  THE  FOLLOWING' ,/, 'LIST: ',//,' 1,  FOR  KEYBOARD  INPUT  VIA  NAME00000192 
>LIST ' , / , ' 2 ,  FOR  FORMATTED  DISK  FILE  INPUT',/, '3,  FOR  INPUT  OF  MINI00000193 
>MUM,  MAXIMUM,  AND  INTERVAL  VALUES ',/,' 4 ,  TO  USE  THE  PREVIOUSLY  PRE00000194 
>SCRIBED  X  VALUES')  00000195 

FORMAT ( II )  00000196 

FORMAT (/,' INPUT  THE  DATA  SET  X  VALUES  (',13,'  MAX)  BY  NAMELIST  WHE00000197 
>RE :',//,' X  =  ARRAY  OF  X  VALUES ',//,' CURRENT  VALUES  ARE:')  00000198 

FORMAT (/, '$PARM  X= _ , _ , _ $',/)  00000199 

FORMAT (/,' DATA  SETS  ARE  INPUT  FROM  FORMATTED  DISK  FILES  AS  X  VALUE00000200 
>S  WHERE: ',//, 'X  =  X  VALUE')  00000201 

FORMAT (/,' INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE  DATA  SET. 00000202 
>')  00000203 

FORMAT (A80)  00000204 

FORMAT (/,' INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES)')  00000205 

FORMAT (' NOTE :  AN  EXAMPLE  FORMAT  IS  " (E15 . 6 )"',/, T7 ,' ENTER  "(*)"  F000000206 
>R  A  FREE  FIELD  READ  (DEFAULT  FORMAT)')  00000207 

FORMAT (/,' INPUT  THE  X  DATA  MAXIMUM,  MINIMUM,  AND  INTERVAL  VALUES  B00000208 
>Y  NAMELIST  WHERE :',//,' XINT  =  X  DATA  POINT  INTERVAL ',/,' XMAX  =00000209 
>  MAXIMUM  X  DATA  VALUE ',/,' XMIN  =  MINIMUM  X  DATA  VALUE ',//,' CURRE000002 10 
>NT  VALUES  ARE:')  00000211 

FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  PROGRAM  LSPRWC')  00000212 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  OPEN  ERROR  ON  UNIT  4,  STATUS  =  00000213 
>"OLD" ,  FILE  =',/)  00000214 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  READ  ERROR  ON  UNIT  4,  FILE  =',/00000215 
>)  00000216 

FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  CLOSE  ERROR  ON  UNIT  4,  FILE  =',00000217 


>/) 

FORMAT ( 'SUBROUTINE  NDWUS  WAS  CALLED  FROM  PROGRAM  LSPRWC') 
FORMAT ( / , ' ERROR  IN  PROGRAM  LSPRWC :  UNFORMATTED  WRITE  ERROR 
>  2') 

FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  PROGRAM  LSPRWC') 
FORMAT (/, T19, '-  ENTER/RETURN  TO  CONTINUE  -') 

FORMAT (Al) 

END 


00000218 
00000219 
ON  UNIT00000220 
00000221 
00000222 
00000223 
00000224 
00000225 
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c 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c* 

c* 

c* 

c* 

c* 

c 

c 


SUBROUTINE  XYDIS(*,*) 


X-Y  DATA  INPUT  SUBROUTINE  (XYDIS) 
REVISION  DATE:  1  JULY  2011 


********************************************************************** 

* 
* 
* 
* 
* 

********************************************************************** 

********************************************************************** 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

********************************************************************** 


SUBROUTINE  XYDIS  DEFINES  THE  X,Y  DATA  POINTS  TO  BE  FIT  FOR  PROGRAM 
LSPRWC . 

INPUT/OUTPUT  VARIABLES: 


X 

Y 


=  ARRAY  OF  X  COORDINATES 
=  ARRAY  OF  Y  COORDINATES 


OUTPUT  VARIABLES: 

NDATA  =  NUMBER  OF  X,Y  DATA  POINTS  TO  BE  FIT 
PARAMETERS : 

MNXYDP  =  MAXIMUM  NUMBER  OF  (X,Y)  DATA  PAIRS  FOR  NAMELIST  INPUT 


IMPLICIT  REAL*8 (A-H,0-Z) 

CH  ARACTE  RCCV*l,Crdimg*80,DFmt*3,File*80,Fmt*80,Fmt_s*80 

PARAMETER  (R_NaN=-999 .E+00) 

PARAMETER  (MNXYDP=101) 

COMMON / NDATA/ NDATA 
COMMON/XINT/XINT 
COMMON / XMAX / XMAX 
COMMON/XMIN/XMIN 

DIMENSION  X (MNXYDP) ,Y (MNXYDP) 

DATA  DFmt/ ' (*) '/ 

NAMELIST/PARM/X, Y 


********************************************************************** 

* 

INPUT  THE  X-Y  DATA  POINTS  TO  BE  FITTED  * 

* 

********************************************************************** 


WRITE (UNIT=6 , FMT=203 ) 
READ (UNIT=5 , FMT=204 ) LCV1 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 
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REWIND  1 

IF (LCV1 . EQ. 2 )  GO  TO  106 


00000061 
00000062 

C  00000063 

C  00000064 

C***********************************************************************00000065 

C*  *00000066 

C*  NAMELIST  INPUT  *00000067 

C*  *00000068 


C 

C 

DO  101  1=1 ,MNXYDP 
X ( I ) =R_NaN 

101  Y ( I ) =0 . 0 

102  WRITE (UNIT=6 , FMT=205 )MNXYDP 
WRITE (UNIT=6 , FMT=206 ) 

CALL  NDRUS ( ' PARM' ,4,*102,*119) 
READ (UNIT=4 , NML=PARM, ERR=102 ) 
CALL  NDRUSE(*102, *119) 

DO  103  NDATA=MNXYDP ,1,-1 
IF(X(NDATA) .NE.R_NaN)  GO  TO  104 

103  CONTINUE 

104  DO  105  1=1 , NDATA 

105  WRITE (UNIT=1 , ERR=120) X ( I ) , Y ( I ) 
GO  TO  112 

C 

C 


00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 


C***********************************************************************00000088 


C*  *00000089 
C*  FORMATTED  DISK  FILE  INPUT  *00000090 
C*  *00000091 

C***********************************************************************00000092 


c 

c 

106 


107 


C 

108 


C 

109 

110 


WRITE (UNIT=6 , FMT=207 ) 

WRITE (UNIT=6 , FMT=208 ) 

READ (UNIT=5 , FMT=209 ) File 
WRITE (UNIT=6 , FMT=210) 

WRITE ( UNIT=6 , FMT=2 1 1 ) 

CALL  DFCUS ( Fmt , DFmt , * 107 ) 

OPEN (UNIT=4 , FILE=File, STATUS= ' OLD ' , ERR=121 ) 

REWIND  4 
NDATA=0 

READ (UNIT=4 , FMT=209 , END=122 , ERR=123 ) Crdimg 
CALL  CBUS01 ( Crdimg, L) 

IF(L.EQ.O)  GO  TO  108 

IF (Fmt (1:3) .NE.DFmt)  READ (Crdimg, FMT=Fmt , ERR=108 )X(1) ,Y(1) 
IF (Fmt (1:3) .EQ.DFmt)  READ( Crdimg, *,ERR=108)X(1) ,Y(1) 
BACKSPACE  4 

IF(Fmt(l:3) .EQ.DFmt)  GO  TO  110 

READ (UNIT=4 , FMT=Fmt , END=111 , ERR=123 ) X ( 1 ) , Y ( 1 ) 

WRITE (UNIT=1 , ERR=120) X ( 1 ) , Y ( 1 ) 

ND  AT  A=ND  AT  A+ 1 
GO  TO  109 

READ (UNIT=4 , FMT=* , END=111 , ERR=123 ) X ( 1 )  ,  Y  ( 1 ) 

WRITE (UNIT=1 , ERR=120) X ( 1 )  ,  Y ( 1 ) 

ND  AT  A=ND  AT  A+ 1 


00000093 

00000094 

00000095 

00000096 

00000097 

00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

00000119 

00000120 
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GO  TO  110 

00000121 

111 

CLOSE (UNIT=4 , STATUS= ' KEEP ' , ERR=125 ) 

00000122 

112 

END  FILE  1 

00000123 

C 

00000124 

C 

00000125 

C***********************************************************************00000126 

c* 

*00000127 

c* 

SEARCH  FOR  THE  MINIMUM  AND  MAXIMUM  X  VALUES  WITHIN  THE  DATA  SET 

*00000128 

c* 

*00000129 

C***********************************************************************00000130 

c 

00000131 

c 

00000132 

REWIND  1 

00000133 

XMAX=-1 . OD+7  5 

00000134 

XMIN=+1 . OD+7  5 

00000135 

XSAV=-1 . OD+7  5 

00000136 

113 

READ (UNIT=1 , END=114 , ERR=124 ) X(l)  ,Y(1) 

00000137 

IF(X(1) . LT.XSAV)  GO  TO  126 

00000138 

XMAX=DMAX1 ( XMAX , X ( 1 ) ) 

00000139 

XMIN=DMIN1 ( XMIN , X ( 1 ) ) 

00000140 

XSAV=X ( 1 ) 

00000141 

GO  TO  113 

00000142 

114 

IF ( XMAX. EQ. XMIN)  GO  TO  127 

00000143 

XINT= (XMAX-XMIN) /I . 0D+02 

00000144 

C 

00000145 

C 

00000146 

C***********************************************************************00000147 

c* 

*00000148 

c* 

DISPLAY  THE  X-Y  DATA  SET 

*00000149 

c* 

*00000150 

C***********************************************************************00000151 

c 

00000152 

c 

00000153 

115 

WRITE (UNIT=6 , FMT=212 ) 

00000154 

CALL  YNOUS (*116, *118, *115) 

00000155 

116 

REWIND  1 

00000156 

WRITE (UNIT=6 , FMT=213 ) 

00000157 

DO  117  1=1 , NDATA 

00000158 

READ (UNIT=1 , ERR=124 ) X(l) ,Y(1) 

00000159 

IF (MOD (I, 20) .NE.O)  GO  TO  117 

00000160 

WRITE (UNIT=6 , FMT=201 ) 

00000161 

READ (UNIT=5 , FMT=202 ) CCV 

00000162 

WRITE (UNIT=6 , FMT=213 ) 

00000163 

117 

WRITE (UNIT=6 , FMT=214 ) I , X ( 1 ) , Y ( 1 ) 

00000164 

WRITE (UNIT=6 , FMT=201 ) 

00000165 

READ (UNIT=5 , FMT=202 ) CCV 

00000166 

118 

RETURN 

00000167 

C 

00000168 

c 

00000169 

£***********************************************************************00000170 

c* 

*00000171 

c* 

ERROR  MESSAGES 

*00000172 

c* 

*00000173 

C***********************************************************************00000174 

c 

00000175 

c 

00000176 

119 

WRITE (UNIT=6 , FMT=215 ) 

00000177 

GO  TO  129 

00000178 

120 

WRITE (UNIT=6 , FMT=216 ) 

00000179 

RETURN  1 

00000180 

B-40 


121 

WRITE (UNIT=6 , FMT=217 ) 

00000181 

CALL  CBUS06 (6, File, *128) 

00000182 

RETURN  1 

00000183 

122 

WRITE (UNIT=6 , FMT=218) 

00000184 

CALL  CBUS06 (6, File, *128) 

00000185 

RETURN  1 

00000186 

123 

WRITE (UNIT=6 , FMT=219 ) 

00000187 

CALL  CBUS06 (6 , File, *128) 

00000188 

RETURN  1 

00000189 

124 

WRITE (UNIT=6 , FMT=220 ) 

00000190 

RETURN  1 

00000191 

125 

WRITE (UNIT=6 , FMT=221 ) 

00000192 

CALL  CBUS06 (6, File, *128) 

00000193 

RETURN  1 

00000194 

126 

WRITE (UNIT=6 , FMT=222 ) 

00000195 

GO  TO  129 

00000196 

127 

WRITE (UNIT=6 , FMT=223 ) 

00000197 

GO  TO  129 

00000198 

128 

WRITE (UNIT=6 , FMT=224 ) 

00000199 

RETURN  1 

00000200 

129 

WRITE (UNIT=6 , FMT=201 ) 

00000201 

READ (UNIT=5 , FMT=202 ) CCV 

00000202 

RETURN  2 

00000203 

C 

00000204 

C 

00000205 

q*  **********************************************************************  qqqqq2Q6 

c* 

*00000207 

c* 

FORMAT  STATEMENTS 

*00000208 

c* 

*00000209 

C***********************************************************************00000210 

c 

00000211 

c 

00000212 

201 

FORMAT (/, T19, ENTER/RETURN  TO  CONTINUE  -') 

00000213 

202 

FORMAT (Al) 

00000214 

203 

FORMAT (/,' SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X-Y  DATA  TO  BE  FITTED00000215 

>  FROM' ,/, 'THE  FOLLOWING  LIST 1 ,  FOR  KEYBOARD  INPUT  VIA 

NAMEL00000216 

>IST ' , / , ' 2 ,  FOR  FORMATTED  DISK  FILE  INPUT NOTE :  THE  X-Y  DATA  M00000217 

>UST  BE  MONATONICALLY  INCREASING  IN  X.') 

00000218 

204 

FORMAT (11) 

00000219 

205 

FORMAT (/,' INPUT  THE  DATA  SET  (X,Y)  PAIRS  (',13,'  MAX)  BY  NAMELIST  00000220 

>WHERE X  =  ARRAY  OF  X  VALUES',/, 'Y  =  ARRAY  OF 

Y  VAL00000221 

>UES ',//,' CURRENT  VALUES  ARE : ' ) 

00000222 

206 

FORMAT (/, '$PARM  X=  ,  ,  , Y=  ,  ,  $',/) 

00000223 

207 

FORMAT (/,' DATA  SETS  ARE  INPUT  FROM  FORMATTED  DISK  FILES  AS  (X 

, Y)  P00000224 

>AIRS  WHERE: ',//, 'X  =  X  VALUE',/, 'Y  =  Y  VALUE') 

00000225 

208 

FORMAT (/,' INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE  DATA 

SET. 00000226 

>') 

00000227 

209 

FORMAT (A80) 

00000228 

210 

FORMAT (/,' INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES)') 

00000229 

211 

FORMAT (' NOTE :  AN  EXAMPLE  FORMAT  IS  " (2E15 . 6 )"',/, T7 ,' ENTER" ( * 

)"  F000000230 

>R  A  FREE  FIELD  READ  (DEFAULT  FORMAT)') 

00000231 

212 

FORMAT (/,' SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION? 

( Y/N) 00000232 

>') 

00000233 

213 

FORMAT (/,T10, 'NO. ' ,T26, 'X' ,T49, 'Y' ) 

00000234 

214 

FORMAT (T8 , 15 , 1P2D23 . 13 ) 

00000235 

215 

FORMAT ( 'SUBROUTINE  NDRUS  WAS  CALLED  FROM  SUBROUTINE  XYDIS') 

00000236 

216 

FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  UNFORMATTED  WRITE  ERROR 

ON  UN00000237 

>IT  1') 

00000238 

217 

FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  OPEN  ERROR  ON  UNIT  4,  STATUS  00000239 

>=  "OLD",  FILE  =' ,/) 

00000240 
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218  FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  UNANTICPATED  EOF  ON  UNIT  4,  00000241 

>FILE  =',/)  00000242 

219  FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  READ  ERROR  ON  UNIT  4,  FILE  ='00000243 

>,/)  00000244 

220  FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  UNFORMATTED  READ  ERROR  ON  UNI00000245 

>T  1')  00000246 

221  FORMAT (/,' ERROR  IN  PROGRAM  LSPRWC:  CLOSE  ERROR  ON  UNIT  4,  FILE  =',00000247 

>/)  00000248 

222  FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  ALL  X-Y  DATA  POINTS  MUST  BE', 00000249 

>/, 'MONOTONICALLY  INCREASING  IN  X.')  00000250 

223  FORMAT (/,' ERROR  IN  SUBROUTINE  XYDIS:  ALL  X-Y  DATA  POINTS  HAVE  THE  00000251 

>SAME  X  VALUE  LEADING  TO',/, 'A  SINGULAR  MATRIX.')  00000252 

224  FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  SUBROUTINE  XYDIS')  00000253 

END  00000254 
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APPENDIX  C 

FORTRAN  TOOL  LIBRARY 


c 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c* 


SUBROUTINE  CBUS01 (Buf ,L)  00000001 

00000002 

00000003 

**********************************************************************00000004 

*00000005 

CHARACTER  BUFFER  UTILITY  SUBROUTINE  (CBUS01)  *00000006 

*00000007 

REVISION  DATE:  3  FEBRUARY  1997  *00000008 

*00000009 

**********************************************************************00000010 

00000011 

**********************************************************************00000012 


c* 

C*  SUBROUTINE  CBUS01  FINDS  THE  NON-TRAILING-BLANK  LENGTH  OF  CHARACTER 
C*  BUFFERS. 

C* 

C*  INPUT/OUTPUT  VARIABLES: 

c* 

C*  Buf  =  CHARACTER  BUFFER 

c* 

C*  OUTPUT  VARIABLES: 

c* 

C*  L  =  NON-TRAILING-BLANK  LENGTH  OF  Buf 

c* 


*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 


C***********************************************************************00000025 


c 

00000026 

c 

00000027 

CHARACTER  Blank*l , Buf * ( * ) 

00000028 

c 

00000029 

SAVE  Blank 

00000030 

c 

00000031 

DATA  Blank/'  '/ 

00000032 

c 

00000033 

DO  101  L=LEN(Buf) ,1,-1 

00000034 

IF(Buf (L:L) .NE. Blank)  RETURN 

00000035 

101 

CONTINUE 

00000036 

L=0 

00000037 

END 

00000038 

C-l 


SUBROUTINE  CBUS06 (LUnit , Buf , * )  00000001 
C  00000002 
C  00000003 
C***********************************************************************00000004 
C*  *00000005 
C*  CHARACTER  BUFFER  UTILITY  SUBROUTINE  (CBUS06)  *00000006 
C*  *00000007 
C*  REVISION  DATE:  11  MARCH  1999  *00000008 
C*  *00000009 
C***********************************************************************00000010 
C  00000011 
C***********************************************************************00000012 
C*  *00000013 
C*  SUBROUTINE  CBUS06  WRITES  CHARACTER  BUFFERS,  LESS  THE  TRAILING  *00000014 
C*  BLANKS,  TO  UNIT  LUnit  USING  DYNAMIC  FORMAT  CONSTRUCTION.  *00000015 
C*  *00000016 
C*  INPUT/OUTPUT  VARIABLES:  *00000017 
C*  *00000018 
C*  LUnit  =  LOGICAL  UNIT  NUMBER  *00000019 
C*  Buf  =  CHARACTER  BUFFER  *00000020 
C*  *00000021 
C***********************************************************************00000022 


c 

c 

c 

c 


CHARACTER  Buf * ( * ) , Fmt*15 

INTEGER  LUnit 

CALL  CBUS01 (Buf ,L) 

IF(L.EQ.O)  L=1 

WRITE ( Fmt , FMT=201 , ERR=101 ) L 

WRITE (UNIT=LUnit , FMT=Fmt , ERR=102 ) (Buf (1:L) ) 

RETURN 


C 

C 


00000023 
00000024 
00000025 
00000026 
00000027 
00000028 
00000029 
00000030 
00000031 
00000032 
00000033 
00000034 
00000035 

C***********************************************************************00000036 

C*  *00000037 

C*  ERROR  MESSAGES  *00000038 

C*  *00000039 

C***********************************************************************00000040 
C  00000041 

C  00000042 

101  WRITE (UNIT=6 , FMT=202 )  00000043 

RETURN  1  00000044 

102  WRITE (UNIT=6 , FMT=2 03 ) LUnit  00000045 

RETURN  1  00000046 

C  00000047 

C  00000048 

C***********************************************************************00000049 
C*  *00000050 

C*  FORMAT  STATEMENTS  *00000051 

C*  *00000052 

C***********************************************************************00000053 
C  00000054 

C  00000055 

201  FORMAT (2H (A, 15 , 1H) )  00000056 

FORMAT (/,' ERROR  IN  SUBROUTINE  CBUS06:  WRITE  ERROR' ,/, 'WHILE  ATTEMP00000057 
>TING  TO  WRITE  VARIABLE  L  TO  Fmt ' ) 


202 

203 


00000058 


FORMAT (/,' ERROR  IN  SUBROUTINE  CBUS06:  WRITE  ERROR' ,/, 'WHILE  ATTEMP00000059 
>TING  TO  WRITE  VARIABLE  Buf  TO  UNIT  =  ',12)  00000060 


C-2 


END 


00000061 


C-3 


SUBROUTINE  DFCUS (Fmt , DFmt , * ) 

C 

C 

Q*********************************************************************** 

c*  * 

C*  DYNAMIC  FORMAT  CONSTRUCTION  UTILITY  SUBROUTINE  (DFCUS)  * 

C*  * 

C*  REVISION  DATE:  20  NOVEMBER  2000  * 

C*  * 

q*********************************************************************** 

c 

Q*********************************************************************** 

c*  * 

C*  SUBROUTINE  DFCUS  READS  A  FORTRAN  FORMAT  FROM  UNIT  5,  STORES  THE  * 

C*  FORMAT  IN  A  CHARACTER  BUFFER,  AND  CHECKS  THE  FORMAT  FOR  ERRORS.  * 

C*  Blank  FORMATS  ARE  REPLACED  WITH  A  DEFAULT  FORMAT.  * 

C*  * 

C*  INPUT  VARIABLES:  * 

C*  * 

C*  DFmt  =  DEFAULT  FORMAT  CHARACTER  BUFFER  * 

C*  * 

C*  OUTPUT  VARIABLES:  * 

C*  * 

C*  Fmt  =  FORMAT  CHARACTER  BUFFER  * 

C*  * 


c 

c 


CHARACTER  CCV*1 ,DFmt* ( * ) , Fmt* ( * ) ,LParn*l , RParn*l 
C 

SAVE  LParn,RParn 
C 

DATA  LParn/ ' ( ' / , RParn/ ' ) ' / 

C 

READ (UNIT=5 , FMT=201 ) Fmt 
ENTRY  DFCUSE( Fmt, DFmt,*) 

CALL  CBUS03 (Fmt) 

CALL  CBUS  0 1 ( Fmt , LFmt ) 

IF(LFmt.GT.O)  GO  TO  101 
CALL  CBUS  0 1 ( DFmt , LDFmt ) 

LFmt=LEN ( Fmt ) 

IF ( LDFmt. GT. LFmt)  GO  TO  102 
Fmt ( 1 : LDFmt ) =DFmt ( 1 : LDFmt ) 


RETURN 

101  IF(Fmt (1 : 1) 

C 

C 

Q**************** 

c* 

c* 

c* 

c 

c 


.EQ. LParn. AND. Fmt (LFmt: LFmt) .EQ. RParn)  RETURN 

******************************************************* 

* 

ERROR  MESSAGES  * 

* 

******************************************************* 


WRITE (UNIT=6 , FMT=202 ) 

CALL  CBUS06 ( 6 , Fmt , *103 ) 

WRITE (UNIT=6 , FMT=203 ) 

READ (UNI T= 5 , FMT=201 ) CCV 
RETURN  1 

102  WRITE ( UNIT=6 , FMT=204 ) LDFmt , LFmt 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 


C-4 


STOP 

103  WRITE (UNIT=6 , FMT=205 ) 
RETURN  1 
C 
C 


00000061 

00000062 

00000063 

00000064 

00000065 


C***********************************************************************00000066 

C*  *00000067 

C*  FORMAT  STATEMENTS  *00000068 

C*  *00000069 

C***********************************************************************00000070 

C  00000071 

C  00000072 

201  FORMAT (A80 )  00000073 

202  FORMAT (/,' ERROR  IN  SUBROUTINE  DFCUS:  ILLEGAL  VALUE  FOR  Fmt  =')  00000074 

203  FORMAT (' FORMATS  MUST  BEGIN  WITH  "("  AND  END  WITH  ,T20 -  EN00000075 

>TER/RETURN  TO  CONTINUE  -')  00000076 

204  FORMAT (/,' FATAL  ERROR  IN  SUBROUTINE  DFCUS THE  NON-TRAILING-B00000077 

>LANK  LENGTH  OF  THE  DEFAULT  FORMAT  CHARACTER' BUFFER,  LDFmt  =',100000078 
>3,'  IS  GREATER  THAN  THE  LENGTH  OF  THE  FORMAT ',/,' CHARACTER  BUFFER, 00000079 
>  LFmt  =' ,13,//, 'EDIT  AND  RECOMPILE  THE  FORTRAN  CODE.')  00000080 

205  FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  SUBROUTINE  DFCUS')  00000081 

END  00000082 


C-5 


c 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c* 

c* 

c* 

c* 

c* 

c 

c 


c 


SUBROUTINE  NDRUS (NBuf , KUnit , * , * ) 


********************************************************************** 

* 

NAMELIST-DIRECTED  READ  UTILITY  SUBROUTINE  (NDRUS)  * 

* 

REVISION  DATE:  23  MARCH  1999  * 

* 

********************************************************************** 

********************************************************************** 

* 

SUBROUTINE  NDRUS  EXECUTES  AN  INTELLIGENT  NAMELIST-DIRECTED  READ  * 

PROCESS  FOR  PROGRAM  PSCCPP  WITH  INPUT  VIA  THE  TEMPORARY  SCRATCH  * 
FILE  ON  LOGICAL  UNIT  NUMBER  Unit.  * 

* 

PARAMETERS :  * 

* 

IBS  =  INPUT  CHARACTER  BUFFER  SIZE  * 

* 

VARIABLES :  * 

* 

NBuf  =  NAMELIST  GROUP-NAME  CHARACTER  BUFFER  * 

IBuf  =  INPUT  CHARACTER  BUFFER  * 

LUnit  =  INPUT  TEMPORARY  SCRATCH  FILE  LOGICAL  UNIT  NUMBER  * 

* 

********************************************************************** 


CHARACTER  AEnd*4 , Amsand*l , Blank*l ,Comma*l ,DEnd*4 ,Dollar*l , IBuf *80 , 
>NBuf * ( * ) , RBuf *80 , Switch*l 

INTEGER  EPontr,SPontr,TPontr 

LOGICAL  OP 

PARAMETER  (IBS=80) 

EQUIVALENCE  (Amsand,AEnd(l : 1) ) , (Dollar, DEnd(l : 1) ) 

SAVE  AEnd , Blank , Comma , DEnd , EPontr , LUnit , SPontr , Switch , TPontr 
DATA  AEnd/ 'SEND'/, Blank/'  '/ , Comma/ ','/, DEnd/ ' $END '/, Switch/ 'U ' / 


********************************************************************** 

* 

CHECK  FOR  A  VALID  NAMELIST  GROUP-NAME  CHARACTER  BUFFER  * 

* 

********************************************************************** 


LUnit=KUnit 

CALL  CBUS04( NBuf, Switch, *110) 
CALL  CBUS03 (NBuf) 

CALL  CBUS01 (NBuf ,L) 

IF(L.EQ.O)  GO  TO  111 
IF (L+2 .GT . IBS )  GO  TO  112 


00000001 

00000002 

00000003 

00000004 

00000005 

00000006 

00000007 

00000008 

00000009 

00000010 

00000011 

00000012 

00000013 

00000014 

00000015 

00000016 

00000017 

00000018 

00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 


C-6 


o  o 


0*********************************************************************** 

c*  * 

C*  OPEN  THE  TEMPORARY  SCRATCH  FILE  ON  LOGICAL  UNIT  NUMBER  LUnit  * 

C*  * 

q*********************************************************************** 

c 

c 

INQUIRE ( UNIT=LUnit , ERR=1 13 , OPENED=OP ) 

IF (OP . EQV. .TRUE. )  CLOSE (UNIT=LUnit , STATUS= ' DELETE' ,ERR=114) 

OPEN (UNIT=LUnit , STATUS= ' SCRATCH ' , FORM= ' FORMATTED ' , ERR=115 ) 

C 

C 

Q*********************************************************************** 

c*  * 

C*  SUPPLY  THE  NAMELIST  INITIALIZATION  AND  NAMELIST  GROUP-NAME  * 

C*  * 

q*********************************************************************** 

c 

c 

CALL  CBUS07 (RBuf ) 

RBuf (2:2) =Dollar 
RBuf ( 3 :L+3 ) =NBuf ( 1 :L) 

WRITE (UNIT=LUnit , FMT=201 , ERR=116 ) RBuf 
C 
C 

Q*********************************************************************** 

c*  * 

C*  READ  THE  NAMELIST  INPUT  BUFFER  * 

C*  * 

q*********************************************************************** 

c 

c 

101  READ (UNIT=5 , FMT=201 , ERR=101 ) IBuf 

C 

C 

Q*********************************************************************** 

c*  * 

C*  CHECK  FOR  A  BLANK  INPUT  BUFFER  * 

C*  * 

q*********************************************************************** 

c 

c 

CALL  CBUS01 ( IBuf , J) 

IF(J.NE.O)  GO  TO  102 
C 
C 

Q*********************************************************************** 

c*  * 

C*  SUPPLY  THE  DEFAULT  CONDITION  * 

C*  * 

q*********************************************************************** 

c 

c 

CALL  CBUS07 ( IBuf) 

IBuf (2:2) =Dollar 

WRITE (UNIT=LUnit , FMT=201 , ERR=116 ) IBuf 
GO  TO  109 


00000061 

00000062 

00000063 

00000064 

00000065 

00000066 

00000067 

00000068 

00000069 

00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 

00000087 

00000088 

00000089 

00000090 

00000091 

00000092 

00000093 

00000094 

00000095 

00000096 

00000097 

00000098 

00000099 

00000100 

00000101 

00000102 

00000103 

00000104 

00000105 

00000106 

00000107 

00000108 

00000109 

00000110 

00000111 

00000112 

00000113 

00000114 

00000115 

00000116 

00000117 

00000118 

00000119 

00000120 
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C***********************************************************************00000121 

C*  *00000122 

C*  SEARCH  FOR  THE  NAMELIST  GROUP-NAME  IN  THE  INPUT  BUFFER  *00000123 

C*  *00000124 

C***********************************************************************00000125 


c 

c 

102 


C 

C 


00000126 
00000127 
00000128 
00000129 
00000130 
00000131 
00000132 
00000133 
00000134 

C***********************************************************************00000135 

C*  *00000136 

C*  CHECK  FOR  A  VALID  NAMELIST  INPUT  *00000137 

C*  *00000138 

C***********************************************************************00000139 


SPontr=l 

CALL  CBUS04(IBuf, Switch, *110) 

CALL  CBUS01 ( IBuf , EPontr ) 

TPontr=INDEX( IBuf (SPontr : EPontr ) , NBuf ( 1 :L) ) 
IF(TPontr.EQ.O)  GO  TO  108 


C 

C 


103 
C 

104 


SPontr=TPontr-l 

IF ( SPontr. LT. 1)  GO  TO  104 

DO  103  1=1, SPontr 

IF(IBuf (1:1) .EQ. Blank)  GO  TO  103 

IF(IBuf (1:1) .EQ. Dollar)  GO  TO  103 

IF(IBuf (1:1) .EQ.Amsand)  GO  TO  103 

RETURN  1 

CONTINUE 

SPontr=TPontr+L 

IF(SPontr.GT. EPontr)  GO  TO  105 

IF(IBuf (SPontr:SPontr) .EQ. Blank)  GO  TO  105 

IF(IBuf (SPontr:SPontr) .EQ. Comma)  GO  TO  105 

IF(IBuf (SPontr:SPontr) .EQ. Dollar)  GO  TO  105 

IF(IBuf(SPontr:SPontr) .EQ.Amsand)  GO  TO  105 

RETURN  1 

SPontr=SPontr-l 
DO  106  1=1, SPontr 
IBuf (I:I)=Blank 
GO  TO  108 


00000140 

00000141 

00000142 

00000143 

00000144 

00000145 

00000146 

00000147 

00000148 

00000149 

00000150 

00000151 

00000152 

00000153 

00000154 

00000155 

00000156 

00000157 

00000158 

00000159 

00000160 

00000161 

00000162 

00000163 

00000164 


C 

105 

106 

C 
C 

C***********************************************************************00000165 

C*  *00000166 

C*  CONTINUE  TO  READ  AND  FLUSH  THE  INPUT  BUFFER  *00000167 

C*  *00000168 

C***********************************************************************00000169 
C  00000170 

C  00000171 

107  READ (UNIT=5 , FMT=201 , ERR=107 ) IBuf  00000172 

CALL  CBUS04(IBuf, Switch, *110)  00000173 

CALL  CBUS01(IBuf, EPontr)  00000174 

108  WRITE (UNIT=LUnit, FMT=2 01 , ERR=116 ) IBuf  00000175 

C  00000176 

C  00000177 

q*  ********************************************************************  *  *00000178 
C*  *00000179 

C*  CHECK  FOR  TERMINATION  IN  THE  INPUT  BUFFER  *00000180 


C-8 


c* 

*00000181 

q*  **********************************************************************  qqqqqI$2 

c 

00000183 

c 

00000184 

SPontr=EPontr-3 

00000185 

IF(IBuf (EPontriEPontr) .EQ. Dollar)  GO  TO  109 

00000186 

IF ( IBuf (EPontr : EPontr ) . EQ . Amsand)  GO  TO  109 

00000187 

IF(IBuf (SPontriEPontr) .EQ.DEnd)  GO  TO  109 

00000188 

IF(IBuf (SPontriEPontr) .EQ.AEnd)  GO  TO  109 

00000189 

GO  TO  107 

00000190 

109 

REWIND  LUnit 

00000191 

RETURN 

00000192 

C 

00000193 

C 

00000194 

C***********************************************************************00000195 

c* 

*00000196 

c* 

CLOSE  THE  TEMPORARY  SCRATCH  FILE  ON  LOGICAL  UNIT  NUMBER  UNIT 

*00000197 

c* 

*00000198 

q*  **********************************************************************  qqqqqI^^ 

c 

00000200 

c 

00000201 

ENTRY  NDRUSE ( * , * ) 

00000202 

CLOSE (UNIT=LUnit , STATUS= ' DELETE ' , ERR=114 ) 

00000203 

RETURN 

00000204 

c 

00000205 

c 

00000206 

q*  **********************************************************************  qqq002Q7 

c* 

*00000208 

c* 

ERROR  MESSAGES 

*00000209 

c* 

*00000210 

Q***********************************************************************QQQQQ2H 

c 

00000212 

c 

00000213 

110 

WRITE (UNIT=6 , FMT=202 ) 

00000214 

RETURN  2 

00000215 

111 

WRITE (UNIT=6 , FMT=203 ) 

00000216 

RETURN  2 

00000217 

112 

WRITE (UNIT=6 , FMT=204 ) 

00000218 

CALL  CBUS06(6,NBuf ,*117) 

00000219 

WRITE (UNIT=6 , FMT=205 ) IBS 

00000220 

RETURN  2 

00000221 

113 

WRITE (UNIT=6 , FMT=206 ) LUnit 

00000222 

RETURN  2 

00000223 

114 

WRITE (UNIT=6,FMT=207) LUnit 

00000224 

RETURN  2 

00000225 

115 

WRITE (UNIT=6 , FMT=208 ) LUnit 

00000226 

RETURN  2 

00000227 

116 

WRITE (UNIT=6 , FMT=209 ) LUnit 

00000228 

RETURN  2 

00000229 

117 

WRITE (UNIT=6 , FMT=210 ) 

00000230 

RETURN  2 

00000231 

C 

00000232 

C 

00000233 

C***********************************************************************00000234 

C*  *00000235 

C*  FORMAT  STATEMENTS  *00000236 

C*  *00000237 

C***********************************************************************00000238 

C  00000239 

C  00000240 


C-9 


201  FORMAT (A80 )  00000241 

202  FORMAT ( 'SUBROUTINE  CBUS04  WAS  CALLED  FROM  SUBROUTINE  NDRUS ' )  00000242 

203  FORMAT (/,' ERROR  IN  SUBROUTINE  NDRUS:  EMPTY  NAMELIST  GROUP-NAME  CHA00000243 

>RACTER' ,/, 'BUFFER,  NBuf ' )  00000244 

204  FORMAT (/,' ERROR  IN  SUBROUTINE  NDRUS:  THE  LENGTH  OF  THE  NAMELIST  GR00000245 

>OUP-NAME  NBuf  =')  00000246 

205  FORMAT ('PLUS  2  EXCEEDS  THE  INPUT  BUFFER  SIZE,  IBS  =',I3)  00000247 

206  FORMAT (/, 'ERROR  IN  SUBROUTINE  NDRUS:  INQUIRE  ERROR  ON  UNIT  ',11)  00000248 

207  FORMAT (/,' ERROR  IN  SUBROUTINE  NDRUS:  CLOSE  ERROR  ON  UNIT  ' , II , / , ' S00000249 

>TATUS  =  "DELETE"')  00000250 

208  FORMAT (/,' ERROR  IN  SUBROUTINE  NDRUS:  OPEN  ERROR  ON  UNIT  ',11,',  ST00000251 

>ATUS  =  "SCRATCH"')  00000252 

209  FORMAT (/,' ERROR  IN  SUBROUTINE  NDRUS:  WRITE  ERROR  ON  UNIT  ',11)  00000253 

210  FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  SUBROUTINE  NDRUS')  00000254 

END  00000255 


C-10 


SUBROUTINE  NDWUS (NBuf , KUnit , * )  00000001 
C  00000002 
C  00000003 

C***********************************************************************00000004 

C*  *00000005 
C*  NAMELIST-DIRECTED  WRITE  UTILITY  SUBROUTINE  (NDWUS)  *00000006 
C*  *00000007 
C*  REVISION  DATE:  23  MARCH  1999  *00000008 
C*  *00000009 
C***********************************************************************00000010 
C  00000011 
C***********************************************************************00000012 


c* 


C*  SUBROUTINE  NDWUS  EXECUTES  AN  INTELLIGENT  NAMELIST-DIRECTED  WRITE 
C*  PROCESS  FOR  PROGRAM  PSCCPP  WITH  OUTPUT  TO  THE  HIGH  LEVEL  PLOT 
C*  COMMAND  FILE  ON  LOGICAL  UNIT  NUMBER  LUnit. 

C* 

C*  PARAMETERS: 

C* 

C*  FieldL  =  FIELD  LENGTH 

C*  OBS  =  OUTPUT  CHARACTER  BUFFER  SIZE 

C* 

C*  VARIABLES: 


C* 

C*  AFL 
C*  CBuf 
C*  DBuf 
C*  IBuf 
C*  N 
C*  N1 
C*  N2 
C*  NBuf 
C*  NNFL 
C*  OBuf 
C*  RBuf 
C*  RFL 
C*  LUnit 
C*  VBuf 
C*  VNFL 
C* 


AVAILABLE  FIELD  LENGTH 

CHARACTER  BUFFER 

DOUBLE  PRECISION  ARRAY  BUFFER 

INTEGER  ARRAY  BUFFER 

NUMBER  OF  ELEMENTS  IN  RBuf  OR  IBuf 

STARTING  ELEMENT  IN  RBuf  OR  IBuf 

ENDING  ELEMENT  IN  RBuf  OR  IBuf 

NAMELIST  GROUP-NAME  CHARACTER  BUFFER 

NAMELIST  NAME  FIELD  LENGTH 

OUTPUT  CHARACTER  BUFFER 

REAL  ARRAY  BUFFER 

REQUIRED  FIELD  LENGTH 

OUTPUT  FILE  LOGICAL  UNIT  NUMBER 

VARIABLE  NAME  CHARACTER  BUFFER 

VARIABLE  NAME  FIELD  LENGTH 


C 

C 

CHARACTER  Apos*2 ,CBuf * ( * ) , Dollar* 1 , Equal* 1 , NBuf* ( * ) ,OBuf *75 , 
>VBuf * ( * ) 

C 

REAL* 8  DBuf 


C 

C 

C 

C 


INTEGER  AFL, EPontr, FieldL, OBS, RFL, SPontr , TPontr , VNFL 
PARAMETER  (FieldL=15 ,OBS=75 ) 

DIMENSION  DBuf (N) , IBuf (N) , RBuf (N) 


EQUIVALENCE  (AFL, NNFL, VNFL) 

C 

SAVE  AFL, Apos, Dollar, EPontr, Equal, ICV, LUnit, OBuf , RFL, SPontr, 
>TPontr , VNFL 
C 

DATA  Apos/" Dollar/ '$'/, Equal/ '='/ 


*00000013 

*00000014 

*00000015 

00000016 

*00000017 

*00000018 

*00000019 

*00000020 

*00000021 

*00000022 

*00000023 

*00000024 

*00000025 

*00000026 

*00000027 

*00000028 

*00000029 

*00000030 

*00000031 

*00000032 

*00000033 

*00000034 

*00000035 

*00000036 

*00000037 

*00000038 

*00000039 

*00000040 

*00000041 

00000042 

00000043 

00000044 

00000045 

00000046 

00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

00000053 

00000054 

00000055 

00000056 

00000057 

00000058 

00000059 

00000060 


C-ll 


o  o  o 


C  00000061 

C  00000062 

C*  **********************************************************************  00000063 
C*  *00000064 

C*  OPEN  THE  NAMELIST  OUTPUT  BUFFER  *00000065 

C*  *00000066 

C***********************************************************************00000067 


c 

c 


LUnit=KUnit 
CALL  CBUS03 (NBuf ) 

CALL  CBUS01 (NBuf ,L) 
NNFL=FieldL 

IF( (L+2) . GT.NNFL)  GO  TO  124 
CALL  CBUS07 (OBuf ) 

SPontr=2 

EPontr=2 

OBuf (SPontr:EPontr)=Dollar 

SPontr=EPontr+l 

EPontr=EPontr+L 

OBuf (SPontr : EPontr ) =NBuf ( 1  :L) 

TPontr=FieldL 

CALL  CBUS01 (OBuf ,L) 

VNFL=TPontr-L-l 

RETURN 


C 

C 


00000068 

00000069 

00000070 

00000071 

00000072 

00000073 

00000074 

00000075 

00000076 

00000077 

00000078 

00000079 

00000080 

00000081 

00000082 

00000083 

00000084 

00000085 

00000086 


00000087 

C***********************************************************************00000088 

C*  *00000089 

C*  INSERT  A  CHARACTER  VARIABLE  INTO  THE  OUTPUT  BUFFER  *00000090 

C*  *00000091 

C***********************************************************************00000092 


c 

c 


101 

102 


00000093 
00000094 
00000095 
00000096 
00000097 
00000098 
00000099 
00000100 
00000101 
00000102 
00000103 
00000104 
00000105 
00000106 
00000107 
00000108 
00000109 
00000110 
00000111 
00000112 
00000113 
00000114 
00000115 
00000116 
00000117 
00000118 
00000119 

***********************************************************************00000120 


ENTRY  NDWUSC ( VBuf ,CBuf , * ) 

CALL  CBUS03 (CBuf ) 

CALL  CBUS01 (CBuf ,L) 

RFL=L+4 
ICV=1 
GO  TO  118 

CALL  CBUS01 (CBuf ,L) 

AFL=0 

TPontr=TPontr+FieldL 
AFL=AFL+FieldL 
IF(RFL.GT.AFL)  GO  TO  102 
IF ( TPontr . GT . OBS )  TPontr=OBS 
SPontr=TPontr-L-2 
EPontr=SPontr 

OBuf ( SPontr : EPontr ) =Apos (1:1) 

SPontr=EPontr+l 

EPontr=TPontr-2 

OBuf (SPontr : EPontr )=CBuf(l:L) 

SPontr=TPontr-l 

EPontr=TPontr 

OBuf (SPontr : EPontr )=Apos 

VNFL=0 

RETURN 


C-12 


C*  *00000121 
C*  INSERT  AN  INTEGER  ARRAY  INTO  THE  OUTPUT  BUFFER  *00000122 
C*  *00000123 

C***********************************************************************00000124 


c 

c 

ENTRY  NDWUSI ( VBuf , N , N1 , N2 , IBuf , * ) 

RFL=FieldL 

ICV=2 

GO  TO  118 

103  DO  107  I=N1 ,N2 

104  TPontr=TPontr+FieldL 

IF ( TPontr . LE . OBS )  GO  TO  106 

JCV=1 

GO  TO  123 

105  GO  TO  104 

106  SPontr=TPontr-FieldL+l 
EPontr=TPontr 

107  WRITE (OBuf (SPontr : EPontr ) , FMT=201 ) IBuf ( I ) 
VNFL=0 

RETURN 

C 

C 


00000125 

00000126 

00000127 

00000128 

00000129 

00000130 

00000131 

00000132 

00000133 

00000134 

00000135 

00000136 

00000137 

00000138 

00000139 

00000140 

00000141 

00000142 

00000143 


C***********************************************************************00000144 


C*  *00000145 
C*  INSERT  A  REAL  ARRAY  INTO  THE  OUTPUT  BUFFER  *00000146 
C*  *00000147 


C***********************************************************************00000148 


c 

c 

ENTRY  NDWUSR ( VBuf , N , N1 , N2 , RBuf , * ) 

RFL=FieldL 

ICV=3 

GO  TO  118 

108  DO  112  I=N1 ,N2 

109  TPontr=TPontr+FieldL 

IF ( TPontr . LE . OBS )  GO  TO  111 

JVC=2 

GO  TO  123 

110  GO  TO  109 

111  SPontr=TPontr-FieldL+l 
EPontr=TPontr 

112  WRITE (OBuf (SPontr : EPontr ) , FMT=202 ) RBuf ( I ) 
VNFL=0 

RETURN 

C 

C 


00000149 

00000150 

00000151 

00000152 

00000153 

00000154 

00000155 

00000156 

00000157 

00000158 

00000159 

00000160 

00000161 

00000162 

00000163 

00000164 

00000165 

00000166 

00000167 


C***********************************************************************q0000168 


C*  *00000169 
C*  INSERT  A  DOUBLE  PRECISION  ARRAY  INTO  THE  OUTPUT  BUFFER  *00000170 
C*  *00000171 


C***********************************************************************00000172 


C  00000173 

C  00000174 

ENTRY  NDWUSD ( VBuf ,N,N1 ,N2 , DBuf , * )  00000175 

RFL=FieldL  00000176 

ICV=4  00000177 

GO  TO  118  00000178 

113  DO  117  I=N1 ,N2  00000179 

114  TPontr=TPontr+FieldL  00000180 
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IF ( TPontr . LE . OBS )  GO  TO  116 

00000181 

JVC=3 

00000182 

GO  TO  123 

00000183 

115 

GO  TO  114 

00000184 

116 

SPontr=TPontr-FieldL+l 

00000185 

EPontr=TPontr 

00000186 

117 

WRITE (OBuf (SPontr : EPontr ) , FMT=202 ) DBuf ( I ) 

00000187 

VNFL=0 

00000188 

RETURN 

00000189 

C 

00000190 

C 

00000191 

C***********************************************************************00000192 

c* 

*00000193 

c* 

INSERT  A  VARIABLE  NAME  INTO  THE  OUTPUT  BUFFER 

*00000194 

c* 

*00000195 

C***********************************************************************00000196 

c 

00000197 

c 

00000198 

118 

CALL  CBUS03 (VBuf ) 

00000199 

CALL  CBUS01 (VBuf ,L) 

00000200 

IF( (L+2)+RFL.GT.OBS)  RETURN  1 

00000201 

119 

IF( (L+2) .LE.VNFL)  GO  TO  120 

00000202 

TPontr=TPontr+FieldL 

00000203 

VNFL=VNFL+FieldL 

00000204 

GO  TO  119 

00000205 

120 

IF ( TPontr+RFL . LE . OBS )  GO  TO  122 

00000206 

JVC=4 

00000207 

GO  TO  123 

00000208 

121 

GO  TO  119 

00000209 

122 

SPontr=TPontr-L-l 

00000210 

EPontr=TPontr-2 

00000211 

OBuf (SPontr : EPontr )=VBuf(l:L) 

00000212 

SPontr=EPontr+2 

00000213 

EPontr=TPontr 

00000214 

OBuf (SPontr:EPontr)=Equal 

00000215 

GO  TO  (101,103,108,113) ,ICV 

00000216 

C 

00000217 

C 

00000218 

q*  **********************************************************************  qqqqq219 

c* 

*00000220 

c* 

FLUSH  THE  OUTPUT  BUFFER 

*00000221 

c* 

*00000222 

q*  **********************************************************************  qqqqq223 

c 

00000224 

c 

00000225 

123 

WRITE (UNIT=LUnit , Fmt=203 , ERR=129 ) OBuf 

00000226 

CALL  CBUS07 (OBuf) 

00000227 

TPontr=0 

00000228 

VNFL=0 

00000229 

GO  TO  (105,110,115,121) , JVC 

00000230 

C 

00000231 

C 

00000232 

q*  **********************************************************************  qqqqq233 

c* 

*00000234 

c* 

CLOSE  THE  NAMELIST  OUTPUT  BUFFER 

*00000235 

c* 

*00000236 

C***********************************************************************00000237 

c 

00000238 

c 

00000239 

ENTRY  NDWUSE ( * ) 

00000240 
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CALL  CBUS01 (OBuf ,L)  00000241 

SPontr=L  00000242 

EPontr=L  00000243 

OBuf (SPontr:EPontr)=Dollar  00000244 

WRITE (UNIT=LUnit , Fmt=203 , ERR=129 ) OBuf  00000245 

IF(LUnit.EQ.6)  WRITE (UNIT=LUnit , Fmt=204 , ERR=129 )  00000246 

RETURN  00000247 

C  00000248 

C  00000249 


q*  **********************************************************************  qqqqq25Q 

c* 

*00000251 

c* 

ERROR  MESSAGES 

*00000252 

c* 

*00000253 

C***********************************************************************00000254 

c 

00000255 

c 

00000256 

124 

WRITE (UNIT=6 , FMT=205 ) 

00000257 

CALL  CBUS06 (6 ,NBuf , *130) 

00000258 

WRITE (UNIT=6 , FMT=206 )NNFL 

00000259 

RETURN  1 

00000260 

125 

WRITE (UNIT=6 , FMT=207 ) 

00000261 

CALL  CBUS06(6, VBuf, *130) 

00000262 

WRITE (UNIT=6 , FMT=208 ) 

00000263 

CALL  CBUS06 ( 6 , CBuf , *130) 

00000264 

WRITE (UNIT=6 , FMT=209 )OBS 

00000265 

RETURN  1 

00000266 

126 

WRITE (UNIT=6 , FMT=207 ) 

00000267 

CALL  CBUS06(6, VBuf, *130) 

00000268 

WRITE (UNIT=6 , FMT=210 )OBS 

00000269 

RETURN  1 

00000270 

127 

WRITE (UNIT=6 , FMT=207 ) 

00000271 

CALL  CBUS06(6, VBuf, *130) 

00000272 

WRITE (UNIT=6 , FMT=211 )OBS 

00000273 

RETURN  1 

00000274 

128 

WRITE (UNIT=6 , FMT=207 ) 

00000275 

CALL  CBUS06 (6 , VBuf , *130) 

00000276 

WRITE (UNIT=6,FMT=212)OBS 

00000277 

RETURN  1 

00000278 

129 

WRITE (UNIT=6,FMT=213)LUnit 

00000279 

RETURN  1 

00000280 

130 

WRITE (UNIT=6 , FMT=214 ) 

00000281 

131 

RETURN  1 

00000282 

C 

00000283 

C 

00000284 

C***********************************************************************q0000285 

c* 

*00000286 

c* 

FORMAT  STATEMENTS 

*00000287 

c* 

*00000288 

C***********************************************************************00000289 

C  00000290 

C  00000291 

00000292 
00000293 
00000294 
00000295 

ROUTINE  NDWUS:  THE  LENGTH  OF  THE  NAMELIST  GR00000296 

00000297 

MELIST  NAME  FIELD  LENGTH  =',I3)  00000298 

RROR  IN  SUBROUTINE  NDWUS:  THE  LENGTH  OF  THE  VARIABLE  NA00000299 

>ME  VBuf  =')  00000300 


201 

FORMAT (I 14, ' , ' ) 

202 

FORMAT (E14. 6, ' , ' 

) 

203 

FORMAT (A7 5) 

204 

FORMAT ( ) 

205 

FORMAT ( / , 'ERROR 

IN 

>OUP-NAME  NBuf  =' 

) 

206 

FORMAT ( 'EXCEEDS 

THE 

207 

FORMAT ( / , ' ERROR 

IN 
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208  FORMAT ('PLUS  THE  LENGTH  OF  THE  CHARACTER  VARIABLE  CBuf  =')  00000301 

209  FORMAT ( 'EXCEEDS  THE  OUTPUT  BUFFER  SIZE,  OBS  =',I3)  00000302 

210  FORMAT ('PLUS  THE  FIELD  LENGTH  FOR  AN  INTEGER  VARIABLE  EXCEEDS  THE  00000303 

>OUTPUT  BUFFER' ,/, 'SIZE,  OBS  =',I3)  00000304 

211  FORMAT ( ' PLUS  THE  FIELD  LENGTH  FOR  A  REAL  VARIABLE  EXCEEDS  THE  OUTP00000305 

>UT  BUFFER' ,/, 'SIZE,  OBS  =',I3)  00000306 

212  FORMAT ('PLUS  THE  FIELD  LENGTH  FOR  A  DOUBLE  PRECISION  VARIABLE  EXCE00000307 

>EDS  THE  OUTPUT  BUFFER' ,/,' SIZE,  OBS  =',I3)  00000308 

213  FORMAT (/,' ERROR  IN  SUBROUTINE  NDWUS :  WRITE  ERROR  ON  UNIT  ',11)  00000309 

214  FORMAT ( 'SUBROUTINE  CBUS06  WAS  CALLED  FROM  SUBROUTINE  NDWUS')  00000310 

END  00000311 
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SUBROUTINE  YNOUS (*,*,*) 

C 

C 

Q********************************************************************* 

c* 

C*  YES/NO/OTHER  UTILITY  SUBROUTINE  (YNOUS) 

C* 

C*  REVISION  DATE:  29  SEPTEMBER  2000 

C* 

Q********************************************************************* 

C 

Q********************************************************************* 

c* 

C*  SUBROUTINE  YNOUS  READS  ONE  BYTE  FROM  UNIT  5  ANTICIPATING  A  'Y'  OR 
C*  'y'  YES  RESPONSE  OR  ALTERNATIVELY  A  'N'  OR  'n'  NO  RESPONSE.  A  YES 
C*  RESPONSE  RESULTS  IN  A  RETURN  1;  A  NO  RESPONSE  RESULTS  IN  A 
C*  RETURN  2.  ALL  OTHER  RESPONSES  RESULT  IN  A  RETURN  3. 

C* 

Q********************************************************************* 

c 

c 

CHARACTER  CCV* 1 , No* 1 , Switch* 1 , Yes  * 1 
C 

SAVE  No, Yes 
C 

DATA  No/ ' N ' / , Switch/ ' U ' / , Yes / ' Y ' / 

C 

C 

Q********************************************************************* 

c* 

C*  READ  ONE  BYTE  AND  COMPARE  WITH  Yes /No 

C* 

Q********************************************************************* 

C 

c 

READ (UNI T= 5 , FMT=201 ) CCV 
CALL  CBUS04(CCV, Switch, *101) 

IF (CCV. EQ. Yes)  RETURN  1 
IF (CCV. EQ . No)  RETURN  2 
RETURN  3 
C 
C 

Q********************************************************************* 

c* 

C*  ERROR  MESSAGES 

C* 

Q********************************************************************* 

c 

c 

101  WRITE (UNIT=6 , FMT=202 ) 

RETURN  3 
C 

Q********************************************************************* 

c* 

C*  FORMAT  STATEMENTS 

C* 

Q********************************************************************* 

c 

c 

201  FORMAT (Al) 


00000001 

00000002 

00000003 

**00000004 

*00000005 

*00000006 

*00000007 

*00000008 

*00000009 

**00000010 

00000011 

**00000012 

*00000013 

*00000014 

*00000015 

*00000016 

*00000017 

*00000018 

**00000019 

00000020 

00000021 

00000022 

00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

**00000029 

*00000030 

*00000031 

*00000032 

**00000033 

00000034 

00000035 

00000036 

00000037 

00000038 

00000039 

00000040 

00000041 

00000042 

**00000043 

*00000044 

*00000045 

*00000046 

**00000047 

00000048 

00000049 

00000050 

00000051 

00000052 

**00000053 

*00000054 

*00000055 

*00000056 

**00000057 

00000058 

00000059 

00000060 


C-17 


00000061 

00000062 


202  FORMAT ( 'SUBROUTINE  CBUS04  WAS  CALLED  FROM  SUBROUTINE  YNOUS ' ) 

END 
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APPENDIX  D 
EXAMPLE  1 


In  this  example,  the  LSPRWC  code  was  applied  to  a  set  of  observations  taken  for  an  object 
accelerating  from  rest  at  a  constant  9.81  m/s2,  the  acceleration  of  gravity.  The  true  behavior  for 
the  object,  in  terms  of  displacement  as  a  function  of  time,  is  given  by 

^2r 

=  constant 
dt 2 

(1) 

=  a, 

(2) 

and  integrating  equation  (1)  twice  gives 

f  =  v0  +  at, 

(3) 

and 

=  v(0. 

(4) 

,  ,  1  2 
x(t)  -  x0  +  v0t  +  —  at 

(5) 

where 

a  =  acceleration, 
t  =  time, 
v  =  velocity, 
v  o  =  initial  velocity, 
x  =  displacement, 
x  o  =  initial  position. 

Then  for 

o  o 

II  II 

O  O 

(6) 

(7) 

equations  (3),  (4),  and  (5)  reduce  to 

v(t)  =  at 

(8) 

and 

x(t)  =  ^ at 2 . 

(9) 

D-l 

The  observations  for  this  example  are  given  in  Table  D-l  along  with  the  true  values  for 
displacement,  velocity,  and  acceleration.  Figures  D-l  through  3  present  X  Window  System 
screen  shots  from  LSPRWC  code  execution  using  this  data  set. 


Table  D-l.  Trajectory  Data 


Time 

(s) 

0+000000E+00 

0+100000E+00 

0+200000E+00 

0+300000E+00 

0.400000E+00 

0.500000E+00 

0tG00000E+00 

0.700000E+00 

0.800000E+00 

0t900000E+00 

0.100000E+01 

0.110000E+01 

0*120000E+01 

0+130000E+01 

0+140000E+01 

0+142784E+01 


True 

Position 

(m) 

0.000000E+00 

0.490500E-01 

0.196200E+00 

0.441450E+00 

0.784800E+00 

(U22G25E+01 

(U7G580E+01 

0.240345E+01 

0.313920E+01 

0.397305E+01 

0.490500E+01 

0.593505E+01 

0.706320E+01 

0.828945E+01 

0+961380E+01 

0+100000E+02 


True 

Velocity 

(m/s) 

0+000000E+00 

0+981000E+00 

0+196200E+01 

0+294300E+01 

0+392400E+01 

0+490500E+01 

0.588G00E+01 

0.G8G700E+01 

0.784800E+01 

0+882900E+01 

0.981000E+01 

(U07910E+02 

0.117720E+02 

(U27530E+02 

0+137340E+02 

0+140071E+02 


True 

Acceleration 

(m/s**2) 

0+981000E+01 

0+981000E+01 

0+981000E+01 

0+981000E+01 

0+981000E+01 

0+981000E+01 

0.981000E+01 

0.981000E+01 

0.981000E+01 

0+981000E+01 

0.981000E+01 

0.981000E+01 

0.981000E+01 

0+981000E+01 

0+981000E+01 

0+981000E+01 


Measured 

Position 

(m) 

0+000000E+00 

0+841G78E+00 

0+262G22E+00 

-0+293G70E+00 

0+177G28E+01 

0+179120E+01 

0.223221E+01 

0.211730E+01 

0.3224GGE+01 

0.306116E+01 

0.466451E+01 

0+669913E+01 

0.781203E+01 

0+790575E+01 

0+105721E+02 

0+920335E+01 
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PROGRAM  LSPRWC  IS  AN  INTERACTIVE  FORTRAN  PROGRAM  TO  PERFORM  A  LEAST 
SQUARES  POLYNOMIAL  REGRESSION  WITH  CONTRA I NTS;  THAT  IS,  A  SET  OF  X-Y  DATA 
POINTS  IS  CURVE  FIT  WITH  AN  NP  ORDER  POLYNOMIAL  OF  THE  FORM 

P ( X ) =B0+B1*X+B2*X**2+B3*X**3+ . . . . +BNP*X**NP 

WITH  ANY  POLYNOMIAL  DERIVATIVES,  ZERO  THROUGH  NP,  SPECIFIED  AT  GIVEN  X 
LOCATIONS,  THE  PROCEDURE  FOLLOWED  IS  THE  METHOD  OF  LEAST  SQUARES  USING 
UNDETERMINED  LAGRANGE  MULTIPLIERS. 

AS  AN  INTERACTIVE  PROGRAM,  LSPRWC  IS  SELF-EXPLANATORY  AND  PROMPTS  FOR  THE 
NECESSARY  INFORMATION.  THE  X-Y  DATA  TO  BE  FITTED  MAY  BE  ENTERED  BY  NAMELIST 
OR  READ  FROM  A  FORMATTED  DISC  FILE.  PROGRAM  LSPRWC  WILL  ALSO  EVALUATE  THE 
RESULTANT  LEAST  SQUARES  POLYNOMIAL  AT  PRESCRIBED  VALUES  OF  X  WHICH,  AGAIN, 
MAY  BE  ENTERED  BY  NAMELIST  OR  READ  FROM  A  FORMATTED  DISK  FILE. 

-  ENTER/RETURN  TO  CONTINUE  - 


SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X-Y  DATA  TO  BE  FITTED  FROM 
THE  FOLLOWING  LIST: 

1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 

2,  FOR  FORMATTED  DISK  FILE  INPUT 

NOTE:  THE  X-Y  DATA  MUST  BE  MONATONICALLY  INCREASING  IN  X. 

2 

DATA  SETS  ARE  INPUT  FROM  FORMATTED  DISK  FILES  AS  (X,Y)  PAIRS  WHERE: 

X  =  X  VALUE 

Y  =  Y  VALUE 

INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE  DATA  SET. 
trajectory.dat 

INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES) 

NOTE:  AN  EXAMPLE  FORMAT  IS  "(2E15.6)" 

ENTER" (*)"  FOR  A  FREE  FIELD  READ  (DEFAULT  FORMAT) 
(E15.6,T61,E15.6) 

SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION?  (Y/N) 

y 


NO. 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


X 

O.OOOOOOOOOOOOOE+OO 
1.0000000000000E-01 
2.0000000000000E-01 
3 . 0000000000000E-01 
4.0000000000000E-01 
5.0000000000000E-01 
6 . 0000000000000E-01 
7.0000000000000E-01 
8 . 0000000000000E-01 
9 . 0000000000000E-01 
l.OOOOOOOOOOOOOE+OO 
1.1000000000000E+00 
1 . 2000000000000E+00 
1 . 3000000000000E+00 
1 . 4000000000000E+00 
1 . 4278400000000E+00 


Y 

O.OOOOOOOOOOOOOE+OO 
8.4167800000000E-01 
2 . 6262200000000E-01 
-2 . 93G7000000000E-01 
1 . 7762800000000E+00 
1 . 7912000000000E+00 
2 . 2322100000000E+00 
2 . 1173000000000E+00 
3 . 224GG00000000E+00 
3 . 0611600000000E+00 
4 . 6645100000000E+00 
G . 8991300000000E+00 
7 . 8120300000000E+00 
7 . 9057500000000E+00 
1 . 0572100000000E+01 
9 . 2033500000000E+00 


-  ENTER/RETURN  TO  CONTINUE  - 


Figure  D-l  LSPRWC  Screen  Shot — Part  1 
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ARE  CONSTRAINTS  DESIRED?  (Y/N) 

N 

INPUT  THE  POLYNOMIAL  ORDER  BY  NAMELIST  WHERE: 

NP  =  THE  POLYNOMIAL  ORDER 

NOTE:  NP  MUST  BE  IN  THE  RANGE  OF  0  TO  9. 

THERE  MUST  BE  AT  LEAST  NP+NCONST+1  X-Y  DATA  POINTS* 

THE  MAXIMUM  VALUE  FOR  NP  IS  9. 

CURRENT  VALUES  ARE: 

$PARM  NP  =  3$ 

Np=2$ 

LEAST  SQUARES  POLYNOMIAL 

P(X)=B(0)+B(l)*X+B(2)*X**2+..*.+B(NP)*X**NP 
I  B(I) 

0  2 . 78068850445S2E-01 

1  -3 . 23645G2759090E-01 

2  5 . 0277277787627E+00 

THE  STANDARD  DEVIATION  FOR  THIS  POLYNOMIAL  OF  ORDER  2  IS  7.17016E-01 
-  ENTER/RETURN  TO  CONTINUE  - 


SHOULD  DATA  POINTS  BE  PRESCRIBED  FOR  EVALUATION  OF  THE  LEAST  SQUARES 
POLYNOMIALS?  (Y/N) 

Y 

SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X  DATA  TO  BE  EVALUATED  FROM  THE  FOLLOWING 
LIST: 

1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 

2,  FOR  FORMATTED  DISK  FILE  INPUT 

3,  FOR  INPUT  OF  MINIMUM,  MAXIMUM,  AND  INTERVAL  VALUES 
3 

INPUT  THE  X  DATA  MAXIMUM,  MINIMUM,  AND  INTERVAL  VALUES  BY  NAMELIST  WHERE: 

XINT  =  X  DATA  POINT  INTERVAL 
XMAX  =  MAXIMUM  X  DATA  VALUE 
XMIN  =  MINIMUM  X  DATA  VALUE 

CURRENT  VALUES  ARE: 

$PARM  XINT  =  0.142784E-01,  XMAX  =  0.142784E+01, 

XMIN  =  O.OOOOOOE+OOS 


SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION?  (Y/N) 
n 

SHOULD  THE  RESULTS  OF  THIS  RUN  BE  QUICK-PLOTTED?  (Y/N) 

N 


Figure  D-2.  LSPRWC  Screen  Shot — Part  2 
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ENTER: 

1,  TO  RESTART  THE  PROGRAM 

2,  TO  CHANGE  THE  CONSTRAINTS 

3,  TO  CHANGE  THE  POLYNOMIAL  ORDER 

4,  TO  CHANGE  THE  PRESCRIBED  X  VALUES  FOR  EVALUATION  OF  THE  POLYNOMIAL 

5,  TO  PLOT  THE  RESULTS 

G,  TO  SAVE  THE  RESULTS  ON  A  FORMATTED  DISC  FILE 
7,  TO  STOP 
6 

INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE. 
np2.dat 

SELECT  THE  TYPE  OF  OUTPUT  FOR  THE  FORMATTED  DISK  FILE  FROM  THE  FOLLOWING  LIST: 

1,  FOR  ALL  INPUT  DATA  AND  RESULTS 

2,  FOR  ONLY  THE  EVALUATED  DATA 
2 

INPUT  THE  EVALUATED  DATA  VARIABLE  LIST  FOR  OUTPUT  BY  NAMELIST  WHERE: 

"X"  =  X  VALUES 

"0“  =  Oth  DERIVATIVE  OF  Y  AT  X 

"1"  =  1th  DERIVATIVE  OF  Y  AT  X 

"2"  =  2th  DERIVATIVE  OF  Y  AT  X 

"R"  =  RADIUS  OF  CURVATURE  AT  X 

"S"  =  ARC  LENGTH 

NOTE:  VARIABLES  WILL  BE  LISTED  IN  THE  ORDER  OF  INPUT 
CURRENT  VALUES  ARE: 

$PARM  AList=  "X\"0"$ 

AList="X","0","l","2"$ 

INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES) 

NOTE:  AN  EXAMPLE  FORMAT  IS  "(FIS^^.G)" 

THE  DEFAULT  FORMAT  IS: 

(13E15.6) 


ENTER: 

1,  TO  RESTART  THE  PROGRAM 

2,  TO  CHANGE  THE  CONSTRAINTS 

3,  TO  CHANGE  THE  POLYNOMIAL  ORDER 

4,  TO  CHANGE  THE  PRESCRIBED  X  VALUES  FOR  EVALUATION  OF  THE  POLYNOMIAL 

5,  TO  PLOT  THE  RESULTS 

6,  TO  SAVE  THE  RESULTS  ON  A  FORMATTED  DISC  FILE 

7,  TO  STOP 
7 


Figure  D-3.  LSPRWC  Screen  Shot — Part  3 


Following  along  through  Figure  D-l,  the  user  was  presented  with  an  introduction  to 
LSPRWC  and  then  asked  for  the  data  source.  Here  the  source  was  selected  as  a  formatted  disk 
file  followed  by  the  file  name  and  read  format.  The  (x,  v)  data  pairs — time  and  measured 
position — were  displayed  for  verification  and  match  the  values  given  in  Table  D-l. 


Moving  on  through  Figure  D-2,  the  code  asked  if  problem  constraints  were  desired  and 
none  were  selected.  The  user  was  then  asked  to  enter  the  polynomial  order  which  was  set  to 
np=  2  by  means  of  NAMELIST.  With  the  problem  thus  fully  defined,  program  LSPRWC 
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returned  values  for  the  coefficients  of  the  least  squares  regression  polynomial,  bo,  b  1,  and  bi,  and 
the  standard  deviation  for  the  fit  to  the  ( x ,  v)  data  pairs  as  input.  This  2nd  order  regression 
polynomial  was  then  evaluated  at  specific  values  of  x,  in  this  case  time,  using  the  minimum, 
maximum,  and  interval  option. 

Following  through  Figure  D-3,  menu  options  were  chosen  to  output  the  evaluated  data  to  a 

dv 

formatted  disk  file  writing,  in  order,  each  value  of  x,  y,  and  the  first  and  second  derivatives,  — , 

d2y 

and  —7.  With  that,  the  LSPRWC  code  menu  option  was  selected  for  program  tennination. 

The  results  from  this  particular  LSPRWC  run  are  shown  in  Figure  D-4  as  plots  of 
displacement  and  velocity  as  a  function  of  time  along  with  the  true  values  and  the  measurements 
of  Table  D-l.  The  2nd  order  least  squares  regression  polynomial  fits  quite  well  in  this  case,  as 
might  be  expected  considering  that  the  true  solution  for  displacement,  given  by  equation  (9),  is 
2nd  order  in  time. 
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_ TIME  (s) _ 

Figure  D-4.  Least  Squares  Polynomial  Regression  for  np  =  2 

The  LSPRWC  code  could  have  been  run  with  any  choice  of  polynomial  order  from 
0  through  9  as  was  done  to  produce  the  results  of  Figure  D-5.  This  plot  of  displacement  as  a 
function  of  time  shows  a  reasonable  fit  for  the  regression  polynomials  for  np  =  2  and  3.  Clearly, 
the  representation  begins  to  degrade  with  polynomial  orders  of  4  and  greater. 

Section  V.  LIMITATIONS  notes  that  the  LSPRWC  code  limits  the  polynomial  order  such 
that  0  <  np  <  9.  This  example  serves  as  excellent  justification  for  this  limitation  since  the 
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9th  order  least  squares  regression  polynomial  gives  the  best  fit  to  the  measured  data  in  Figure  D-5 
yet  the  polynomial  representation  is  quite  unlike  the  true  solution. 


().()  0.3  0.6  0.9  1.2  1.5 


0.0  0.3  0.6  0.9  1.2  1.5 

_ TIME  (s) _ 

Figure  D-5.  Least  Squares  Polynomial  Regressions  for  np  =  2,  3,  4,  5,  6,  7,  8,  9 

Examining  once  again  the  plots  of  displacement  and  velocity  as  a  function  of  time  in 
Figure  D-4,  it  is  obvious  that  this  2nd  order  regression  polynomial  does  not  satisfy  either  equation 
(6)  or  (7);  that  is  to  say,  the  displacement  and  velocity  are  non-zero  at  time  equal  zero. 
Furthermore,  the  polynomial  coefficients  bo  and  b\  are  non-zero  as  shown  in  Figure  D-2.  Still, 
the  acceleration  a  =  2  x  b2  =  10.055  m/s2  is  reasonably  close  to  the  true  value  of  9.81  m/s2. 
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Since  the  initial  displacement  and  velocity  are  known,  these  could  be  used  as  constraints 
with  the  LSPRWC  code  to  improve  the  regression.  Figures  D-6  through  9  present^  Window 
System  screen  shots  from  LSPRWC  code  execution  using  again  the  observations  of  Table  D-l 
but  with  these  constraints. 


AS  AN  INTERACTIVE  PROGRAM,  LSPRWC  IS  SELF-EXPLANATORY  AND  PROMPTS  FOR  THE 
NECESSARY  INFORMATION.  THE  X-Y  DATA  TO  BE  FITTED  MAY  BE  ENTERED  BY  NAMELIST 
OR  READ  FROM  A  FORMATTED  DISC  FILE.  PROGRAM  LSPRUC  UILL  ALSO  EVALUATE  THE 
RESULTANT  LEAST  SQUARES  POLYNOMIAL  AT  PRESCRIBED  VALUES  OF  X  UHICH,  AGAIN, 
MAY  BE  ENTERED  BY  NAMELIST  OR  READ  FROM  A  FORMATTED  DISK  FILE. 

-  ENTER/RETURN  TO  CONTINUE  - 


SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X-Y  DATA  TO  BE  FITTED  FROM 
THE  FOLLOWING  LIST: 

1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 

2,  FOR  FORMATTED  DISK  FILE  INPUT 

NOTE:  THE  X-Y  DATA  MUST  BE  MONATONICALLY  INCREASING  IN  X. 

2 

DATA  SETS  ARE  INPUT  FROM  FORMATTED  DISK  FILES  AS  (X,Y)  PAIRS  WHERE: 

X  =  X  VALUE 

Y  =  Y  VALUE 

INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE  DATA  SET. 
trajectory.dat 

INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES) 

NOTE:  AN  EXAMPLE  FORMAT  IS  "(2E15.6)" 

ENTER" (*)"  FOR  A  FREE  FIELD  READ  (DEFAULT  FORMAT) 
(E15.6,T61,E15.6) 

SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION?  (Y/N) 
n 


Figure  D-6.  LSPRWC  Screen  Shot — Part  1 


D-9 


ARE  CONSTRAINTS  DESIRED?  (Y/N) 

y 

SELECT  THE  SOURCE  OF  INPUT  FOR  THE  CONSTRAINTS  FROM  THE  FOLLOWING  LIST: 

1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 

2,  FOR  FORMATTED  DISK  FILE  INPUT 

NOTE:  THE  CONSTRAINTS  ARE  COMPLETELY  INDEPENDENT  OF  THE  X-Y  DATA  POINTS  TO  BE 
FITTED. 

1 

INPUT  THE  CONSTRAINTS  (20  MAX)  BY  NAMELIST  WHERE: 

CV  =  ARRAY  OF  POLYNOMIAL  DERIVATIVE  VALUES 

NC  =  ARRAY  OF  POLYNOMIAL  DERIVATIVE  ORDERS  (CONSTRAINT  ORDER) 

XC  =  ARRAY  OF  POLYNOMIAL  DERIVATIVE  LOCATIONS  (CONSTRAINT  X  VALUES) 

NOTE:  NC  VALUES  MUST  BE  IN  THE  RANGE  0  TO  NP. 

CURRENT  VALUES  ARE: 

$PARM  CV= _ ,NC= _ _ _ XC= _ $ 

XC=2*0.0,NC=0,1,CV=2*0.0$ 

SHOULD  THE  CONSTRAINTS  BE  DISPLAYED  FOR  VERIFICATION?  (Y/N) 

Y 

NO.  XC(I)  NC(I)  CV(I) 

1  O.OOOOOOOOOOOOOD+OO  0  O.OOOOOOOOOOOOOD+OO 

2  O.OOOOOOOOOOOOOD+OO  1  O.OOOOOOOOOOOOOD+OO 

-  ENTER/RETURN  TO  CONTINUE  - 


INPUT  THE  POLYNOMIAL  ORDER  BY  NAMELIST  WHERE: 

NP  =  THE  POLYNOMIAL  ORDER 

NOTE:  NP  MUST  BE  IN  THE  RANGE  OF  0  TO  9. 

THERE  MUST  BE  AT  LEAST  NP+NCONST+1  X-Y  DATA  POINTS. 

THE  MAXIMUM  VALUE  FOR  NP  IS  9. 

CURRENT  VALUES  ARE: 

$PARM  NP  =  3$ 

NP=2$ 

LEAST  SQUARES  POLYNOMIAL 

P(X)=B(0)+B(l)*X+B(2)*X**2+....+B(NP)*X**NP 
I  B(I) 

0  O.OOOOOOOOOOOOOD+OO 

1  O.OOOOOOOOOOOOOD+OO 

2  4 . 96149408B8317D+00 

THE  STANDARD  DEVIATION  FOR  THIS  POLYNOMIAL  OF  ORDER  2  IS  8.03870D-01 
-  ENTER/RETURN  TO  CONTINUE  - 


Figure  D-7.  LSPRWC  Screen  Shot — Part  2 
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SHOULD  DATA  POINTS  BE  PRESCRIBED  FOR  EVALUATION  OF  THE  LEAST  SQUARES 
POLYNOMIALS?  (Y/N) 
a 

SELECT  THE  SOURCE  OF  INPUT  FOR  THE  X  DATA  TO  BE  EVALUATED  FROM  THE  FOLLOWING 
LIST: 

1,  FOR  KEYBOARD  INPUT  VIA  NAMELIST 

2,  FOR  FORMATTED  DISK  FILE  INPUT 

3,  FOR  INPUT  OF  MINIMUM,  MAXIMUM,  AND  INTERVAL  VALUES 
3 

INPUT  THE  X  DATA  MAXIMUM,  MINIMUM,  AND  INTERVAL  VALUES  BY  NAMELIST  WHERE: 

XINT  =  X  DATA  POINT  INTERVAL 
XMAX  =  MAXIMUM  X  DATA  VALUE 
XMIN  =  MINIMUM  X  DATA  VALUE 

CURRENT  VALUES  ARE: 

$PARM  XINT  =  0.142784E-01,  XMAX  =  0.142784E+01, 

XMIN  =  0.000000E+00$ 


SHOULD  THE  X-Y  DATA  BE  DISPLAYED  FOR  VERIFICATION?  (Y/N) 
n 

SHOULD  THE  RESULTS  OF  THIS  RUN  BE  QUICK-PLOTTED?  (Y/N) 
n 

ENTER: 

1,  TO  RESTART  THE  PROGRAM 

2,  TO  CHANGE  THE  CONSTRAINTS 

3,  TO  CHANGE  THE  POLYNOMIAL  ORDER 

4,  TO  CHANGE  THE  PRESCRIBED  X  VALUES  FOR  EVALUATION  OF  THE  POLYNOMIAL 

5,  TO  PLOT  THE  RESULTS 

G,  TO  SAVE  THE  RESULTS  ON  A  FORMATTED  DISC  FILE 
7,  TO  STOP 
6 

INPUT  THE  FILE  NAME  OF  THE  FORMATTED  DISK  FILE. 
np2_c.dat 

SELECT  THE  TYPE  OF  OUTPUT  FOR  THE  FORMATTED  DISK  FILE  FROM  THE  FOLLOWING  LIST: 

1,  FOR  ALL  INPUT  DATA  AND  RESULTS 

2,  FOR  ONLY  THE  EVALUATED  DATA 
2 

INPUT  THE  EVALUATED  DATA  VARIABLE  LIST  FOR  OUTPUT  BY  NAMELIST  WHERE: 

"X"  =  X  VALUES 

"0"  =  Oth  DERIVATIVE  OF  Y  AT  X 

"1"  =  1th  DERIVATIVE  OF  Y  AT  X 

"2"  =  2th  DERIVATIVE  OF  Y  AT  X 

"R"  =  RADIUS  OF  CURVATURE  AT  X 

"S"  =  ARC  LENGTH 

NOTE:  VARIABLES  WILL  BE  LISTED  IN  THE  ORDER  OF  INPUT 
CURRENT  VALUES  ARE: 

$PARM  AList=  "X","0"$ 

AList="X”,"0","l","2"$ 


Figure  D-8.  LSPRWC  Screen  Shot — Part  3 
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INPUT  THE  DATA  FILE  FORMAT  (INCLUDE  PARENTHESES) 

NOTE:  AN  EXAMPLE  FORMAT  IS  "(F15.7,E15.6)" 

THE  DEFAULT  FORMAT  IS: 

(13E15.6) 


ENTER: 

1,  TO  RESTART  THE  PROGRAM 

2,  TO  CHANGE  THE  CONSTRAINTS 

3,  TO  CHANGE  THE  POLYNOMIAL  ORDER 

4,  TO  CHANGE  THE  PRESCRIBED  X  VALUES  FOR  EVALUATION  OF  THE  POLYNOMIAL 

5,  TO  PLOT  THE  RESULTS 

6,  TO  SAVE  THE  RESULTS  ON  A  FORMATTED  DISC  FILE 

7,  TO  STOP 
7 


Figure  D-9.  LSPRWC  Screen  Shot — Part  4 


Following  along  through  Figure  D-6,  the  user  was  presented  with  the  introduction  to 
LSPRWC,  the  data  source  was  selected  as  a  formatted  disk  file,  and  the  file  name  and  read 
format  were  entered. 


Continuing  through  Figure  D-7,  the  code  asked  if  problem  constraints  were  desired  and  a 
positive  response  was  entered  with  the  constraints  defined  by  means  of  NAMELIST.  The 
constraints  were  then  displayed  for  verification  to  confirm  that  the  0th  and  1st  derivatives  were  set 
to  zero  at  time  zero.  The  polynomial  order  was  again  set  to  np  =  2  and  with  the  problem  thus 
fully  defined,  program  LSPRWC  returned  values  for  the  coefficients  of  the  least  squares 
regression  polynomial. 

Following  now  through  Figure  D-8,  this  second-order  regression  polynomial  was  evaluated 
at  specific  values  of  x  again  using  the  minimum,  maximum,  and  interval  option.  Menu  options 
were  then  chosen  to  output  the  evaluated  data  to  a  formatted  disk  file  writing,  in  order,  each 
value  of  x,  y,  and  the  1 st  and  2nd  derivatives. 

Finally,  in  Figure  9,  the  LSPRWC  code  menu  option  was  selected  for  program 
termination. 

The  results  from  this  LSPRWC  run  are  shown  in  Figure  D-10  as  plots  of  displacement  and 
velocity  as  a  function  of  time  along  with  the  true  values  and  the  measurements  of  Table  D-l. 

This  2nd  order  least  squares  regression  polynomial  agrees  exceptionally  well  with  the  true 
relations  for  displacement  and  velocity  having  constraints  imposed  to  satisfy  equations  (6)  and 
(7). 
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TIME  (s) 


Figure  D-10.  Least  Squares  Polynomial  Regression  for  np  =  2  with  Constraints 

Indeed,  the  polynomial  coefficients  bo  and  b\  are  now  zero  valued,  as  shown  in  Figure  D-7. 
Furthermore,  the  acceleration,  a  =  2  x  h2  =  9.92  m/s2,  is  much  closer  to  the  true  value  of 
9.81  m/s  demonstrating  the  value  of  added  constraints  to  the  least  squares  polynomial  regression 
procedure  when  they  are  known. 
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APPENDIX  E 
EXAMPLE  2 


For  this  example,  the  LSPRWC  code  was  applied  to  the  representation  of  a  set  of 
observations  artificially  generated  through  the  application  of  random  dispersions  to  a  hyperbolic 
spiral.  The  hyperbolic  spiral  in  polar  coordinates  is  given  by 

R  =  a/6  (1) 

or  in  Cartesian  coordinates 

x  =  a/cos(6), 
y  =  a/sin(0), 

where 

a  =  constant, 

R  =  radius, 

9  =  angle, 
x  =  ordinate, 
y  =  absicca. 

A  partial  list  of  the  observations  for  this  example  are  given  in  Table  E-l,  along  with  the 
true  values  for  the  x  and  v  coordinates.  Note  that  these  observations  were  generated  using 
equation  (1)  for  equal  1°  increments  of  6  from  90°  through  720°  with  a  dispersion  applied  to  R 
within  a  uniformly  distributed  random  resolution.  These  polar  coordinates  were  then  converted 
to  Cartesian  coordinates  to  produce  the  values  given  in  Table  E-l. 


(2) 

(3) 
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Table  E-l.  Observations,  Cartesian  Coordinates 


True 

X 

0.000000E+00 

-0409885E-01 

-0.217347E-01 

-0.322433E-01 

-0.425186E-01 

-0.525G48E-01 

-0.623858E-01 

-0.719856E-01 

-0.813G77E-01 

-0.905357E-01 

-0.994931E-01 

-0.108243E+00 

-0416789E+00 

-0.125133E+00 

-0.133280E+00 

-0.141231E+00 

-0448989E+00 

-0.156558E+00 

-04G3939E+00 

-0.171135E+00 

-0478148E+00 


0.772813E-01 

0.77G233E-01 

0.77940GE-01 

0.782333E-01 

0.785014E-01 

0.787448E-01 

0.789636E-01 

0.791578E-01 

0.793274E-01 

0.794723E-01 

0.795926E-01 

0.796885E-01 

0.797598E-01 

0.798066E-01 

0.798290E-01 

0.798271E-01 

0.798009E-01 

0.797505E-01 

0.796760E-01 

0.795775E-01 

0.795775E-01 


True 

Y 

0.63GG20E+00 

0.629528E+00 

0.622401E+00 

0.615239E+00 

0.G08045E+00 

0.G00818E+00 

0.593562E+00 

0.586275E+00 

0.5789G1E+00 

0.571G20E+00 

0.564253E+00 

0.55G8G2E+00 

0.549448E+00 

0.542013E+00 

0.534556E+00 

0.527081E+00 

0.519587E+00 

0.512077E+00 

0.504551E+00 

0.497011E+00 

0.489458E+00 


Measured 

X 

0.000000E+00 

-040G533E-01 

-0.207991E-01 

-0.328G93E-01 

-0.434283E-01 

-0.500020E-01 

-0.630342E-01 

-0.722125E-01 

-0.751539E-01 

-0.929371E-01 

-0.984405E-01 

-0.113103E+00 

-0.118323E+00 

-0.128515E+00 

-0.139035E+00 

-0.1431G0E+00 

-0.157430E+00 

-0.1G142GE+00 

-0471451E+00 

-0.169825E+00 

-0.163469E+00 


Measured 

Y 

0.637909E+00 

0.G10330E+00 

0.595609E+00 

0.627184E+00 

0.G21054E+00 

0.571525E+00 

0.599730E+00 

0.588124E+00 

0.534748E+00 

0.586782E+00 

0.558284E+00 

0.5818G5E+00 

0.556665E+00 

0.556658E+00 

0.557639E+00 

0.534281E+00 

0.549024E+00 

0.528002E+00 

0.527G71E+00 

0.493208E+00 

0.449128E+00 


-0467373E-01 

-0.389375E-01 

-0.353019E-01 

-0.323478E-01 

-0413676E-01 

-0.80G209E-02 

-0.292896E-01 

-0.221176E-01 

-0.173756E-01 

-0496431E-01 

-0.129619E-01 

-0.142397E-01 

-0.110859E-01 

-0.494223E-02 

-0.757223E-02 

-0.817499E-02 

-0.220901E-02 

-0.415463E-02 

-0.950917E-03 

0.000000E+00 

0.000000E+00 


-0.266101E-01 
-0.252213E-01 
-0.238288E-01 
-0.224330E-01 
-0.210344E-01 
-0.19G333E-01 
-0.182302E-01 
-0  468255E-01 
-0.154197E-01 
-0440131E-01 
-0.126062E-01 
-0.111995E-01 
-0.979327E-02 
-0.838801E-02 
-0.G98413E-02 
-0.558206E-02 
-0.418219E-02 
-0.278495E-02 
-0439075E-02 
0.000000E+00 
0.000000E+00 


0.48G087E-01 
0419837E+00 
0415467E+00 
0412810E+00 
0.424245E-01 
0.323353E-01 
0  426867E+00 
0404055E+00 
0.893899E-01 
0411402E+00 
0.818384E-01 
0401321E+00 
0.90287GE-01 
0.470222E-01 
0.8G5510E-01 
041G908E+00 
0.421504E-01 
0418973E+00 
0.544780E-01 
0.751481E-01 
0.798177E-01 


The  full  set  of  observations  from  Table  E-l  are  presented  in  the  plot  of  Figure  E-l  and 
show  immediately  that  this  set  of  ( x ,  v)  data  pairs  or  observations  will  not  be  suitable  for  a  least 
squares  polynomial  regression.  The  true  position  data  depicting  the  hyperbolic  spiral,  as 
indicated  in  red,  exhibits  locations  of  infinite  slope  which  cannot  be  represented  with 
polynomials.  Furthermore,  the  (x,y)  observations  may  well  be  multivalued  for  a  given  x 
location,  as  is  certainly  the  case  for  the  hyperbolic  spiral.  Although  artificially  generated 
specifically  for  this  example,  similar  Cartesian  coordinate  data  sets  do  occur  and  pose  similar 
problems  for  regression. 
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Figure  E-l.  Measured  and  True  Position  in  (x,  y)  Coordinates 


Something  must  be  done  to  fit  this  (x,  y)  data  set  of  Table  E-l  with  a  least  squares 
polynomial  regression  and  this  is  accomplished,  at  least  for  this  example,  through  a  coordinate 
transformation.  Figure  E-2  shows  the  data  in  polar  coordinates,  data  which  is  clearly  well 
behaved  and  amenable  to  treatment  with  the  method  of  least  squares  polynomial  regression. 
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Theta  (rad) 

Figure  E-2.  Measured  and  True  Position  in  Polar  Coordinates 

The  LSPRWC  code  was  then  run  with  the  observations  in  the  polar  coordinates  of  Figure 
E-2  for  polynomial  orders  from  2  through  5  to  produce  the  results  of  Figure  E-3.  This  plot  of  R 
versus  6  shows  a  reasonable  fit  for  the  regression  polynomials,  although  none  truly  duplicate  the 
hyperbolic  spiral  a/9  as  shown  in  the  figure.  Perhaps  this  should  not  be  unanticipated  since  R  is  a 
function  of  9  to  the  -1  power,  while  the  LSPRWC  code  limits  the  np  order  regression 
polynomial,  in  this  case,  to 


np 

m  =  Yjb>°J 

7=0 


for  0  <  np  <  9. 
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Figure  E-3.  Least  Squares  Polynomial  Regression  for  np  =  2,  3,  4,  5 


Assuming  the  5  th  order  regression  polynomial  to  give  the  best  fit,  a  set  of  evaluations  from 
that  polynomial  over  the  range  90°  <  6  <  720°  were  converted  to  Cartesian  coordinates  and 
plotted  in  Figure  E-4  along  with  the  observations  and  true  values  from  Table  E-l.  Surprisingly, 
the  lit  is  quite  good  with  the  largest  departure  from  the  theoretical  occurring  for  values  of  9  near 
90°  rather  than  720°,  as  might  be  surmised  from  Figure  E-2. 
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